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ABSTRACT 

Ekman's    linear    equations    for    time  dependent    flow    (neglecting 
wind  stress)    are   solved  using  a    time   dependent   Green's    function   and 
the  method   suggested  by   Welander    (1957).      The   solution   represents 
the  vertical    velocity  profile    in    terms    of   the   local   time  history    of 
the   changes    in   sea    surface   elevation    determined  by    the   divergence 
of   the   flow    in   the   vertically    integrated  continuity   equation. 

A  fully    implicit    finite-difference  scheme    is    developed  to  re- 
present  a    time   dependent    seiche    oscillating   across    a    shallow    infinite 
channel.      The   transients    associated  with    the    formation    of    the  bottom 
spiral   are    clearly   represented  by    the  model   and   the    influence   of 
friction   and   Coriolis   are    individually   and   collectively    introduced. 
The  model   allows    independent    calculation    of   velocity,    volume   trans- 
port,   sea    surface   elevation,    bottom   stress,    and   the   total    energy 
balance   of   the  system.      The  numerical   scheme  provides    a   method   for 
adequately    describing   and    investigating   certain   classes    of   time 
dependent   motion,    and    its    development   suggests   a   mechanism  for 
improving    numerical  prediction    of   storm  surge. 
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I.       INTRODUCTION 


"The  purpose   of   computing    is    insight    not   numbers" 

(Hamming,    1962) 


A.       BACKGROUND 

Oceanographers    use  a   number    of   techniques    for    the   study   of 
oceanic    flow    including   analysis    of   hydrographic    data,    use    of 
hydraulic    laboratory   models,    and  mathematical    treatment    of   the 
basic   equations    of  motion.       Initially    oceanic   circulation   was 
treated  as    a    steady    state  problem,    and   various    investigators    have 
been   successful    in   representing   and   explaining   many    of   the  gross 
features    of    oceanic  circulation   using   each   of   these   techniques. 

Mathematical    treatment    of    oceanic    circulation   began   with 
Ekman's   analysis    (1905)    of  wind-driven   circulation   and  his    results 
have  been    extended  by   a    number    of  authors.       (Sverdrup,    1947,    Munk, 
1950,    Bryan,    1959)         For    the   most  part    the  study    of    large   scale 
ocean    dynamics   has   been   done  with  steady    state   formulations. 

Mathematical   treatment    of   time   dependent   motion    is    far    less 
complete    in   the   literature   than    that    of   steady   state.       It    is   also 
far   more    difficult   to   compare   results    due   to   the  scarcity    of 
reliable   data.      The    complex   nature   of  the   solutions    to   time   de- 
pendent  boundary   value  problems    further    increases    this    difficulty. 
These  solutions    take   the   form   of   varying    differential,    integral, 
and    infinite  series    equations    whose  numerical    evaluation  was    not 
feasible  until    the  advent    of   the  high   speed   computer    and   the   de- 
velopment   of  appropriate   numerical    techniques. 


In  the  last  decade  these  methods  have  been  applied,  through 
numerical  models,  to  the  problems  of  explaining  and  predicting  the 
various  details  of  fluid  motion  in  both  the  atmosphere  and  the 
oceans.   Although  much  success  has  been  attained  with  these  methods, 
many  critical  questions  still  need  explanation.   Many  of  these 
questions  relate  to  the  structure  of  the  various  frictional  boundary 
layers,  the  details  of  their  development,  and  the  resulting  effects 
on  the  interior  flow. 

Many  complicated  engineering  problems  occur  in  near  shore  and 
other  shallow  water  regions.   The  time  dependent  nature  of  these 
areas  is  well  documented  and  the  resulting  flow  is  greatly  modified 
due  to  interaction  with  the  ocean  bottom.   Therefore  an  understanding 
of  the  nature  of  time  dependent  flow  at  all  levels  is  fundamental  to 
dealing  with  problems  in  these  areas.   The  question  of  how  readjust- 
ment takes  place  and  how  the  many  transients  can  be  appropriately 
represented  in  the  bottom  frictional  boundary  layer  has  not  been 
adequately  explained.   Further,  it  is  necessary  to  understand  the 
dynamic  coupling  between  the  flow  and  the  development  of  the  bottom 
boundary  layer  before  accurate  modeling  techniques  can  be  developed. 
It  is  towards  this  understanding  that  this  work  was  directed. 

B.   PREVIOUS  RELATED  RESEARCH 

Work  by  Welander  (1957)  dealt  with  obtaining  time  dependent 
solutions  using  Ekman  dynamics.   His  results  were  presented  in 
terms  of  a  time  dependent  Green's  function  and  were  related  to  the 
problems  of  predicting  tides  and  storm  surges.   He  demonstrated 
that  the  velocity  profile,  volume  transport,  sea  level  elevation, 
and  bottom  stress  could  be  exactly  represented  by  his 


"integr o-dif ferent ial"    solution.       It   was    suggested    that    this 
solution   could  be   used  to    improve  storm   surge   and  tide  predictions. 
A  step-by-step   numerical    scheme  was   proposed,    but    the   details    of 
the  scheme   were  not  presented   in   the   literature. 

Gait    (1970)    covered  much   of   the   same  material    that   Welander 
presented  but    used  a    different   method   of   solution.      The   solution 
allowed  a   more   explicit    formulation    of   results    which  provided  more 
insight    into   the  possible   transients   present  and   the  relevant    time 
scales    of   each.      He   discussed   the   need   for    information   relating    to 
time    dependent    flow    in    shallow  water    regions,    and   suggested   that 
the   numerical    investigation   of   the    interaction   between   a    seiche 
and    its   associated  transient   bottom   spiral   might   be   useful    in   de- 
termining  how   the   various    transients   present   are   effectively 
coupled. 

C.      THESIS   OBJECTIVES 

Since   the  major    area    of   interest   was    to   be   the  bottom  boundary 
layer,    the   first    objective  was    to   formally   solve   the  governing 
equations    neglecting   wind   stress.      The   solution  would   therefore   not 
allow  a   surface  Ekraan    spiral   to   develop   and   effect    the    interior 
flow.      It   will  be  showed  that    the   solution   can   be  represented  by 
two   equations,    one    integral   and   one    differential.      These   equations 
are   similar    in    form   to   the   ones    developed  by    Welander    and  Gait. 
Although   both    of   these   authors   proposed   that    these   equations    could 
be   solved  numerically   using   step-by-step  procedures,    it  was    not 
initially   apparent   whether    a    stable  and  accurate   numerical    scheme 
was    feasible.      Therefore   the   second   objective   was    to   formulate  a 
numerical    model   capable   of   simultaneously    solving    these   equations. 
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Finally,  as  suggested  by  Gait,  this  model  would  be  applied  to 
represent  a  seiche  to  examine  the  relationships  which  developed 
between  the  various  transients. 
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II.       THE   FORMAL    SOLUTION 

A.      GOVERNING   EQUATIONS  AND  ASSUMPTIONS 

The  momentum   equations    considered  are   essentially    those   used 
by    Ekraan.      The   non-linear    and   horizontal    friction    terms    are   assumed 
small   and   neglected.      The  Coriolis   parameter,    the   density,    and   the 
coefficient    of  vertical    friction  are  all   assumed   constant.      Thus: 

2 

bu  -  .     b    u  b£  , ,  . 

fv    _   k_  _    .   g  _-  (1) 

Oz 

__  +   fu   .   k  =   .   g  ^  .  (2) 

Oz 

Where  u  and  v   are   velocities    in   the  x   and   y-dir  ect  ions  ,    f   the 
Coriolis  parameter,    g   the  acceleration    of  gravity,    and  k   the 
constant    coefficient    of   vertical    friction.      ?(x,y,t)    is    the   ele- 
vation  of  the  free  surface   and   the  z-axis    is   considered  positive 
up. 

The   third  equation   necessary    to   solve    for    u,    v,    and  §    is    the 
vertically    integrated   continuity   equation. 


|L  +  |u  +  bV  = 

bt        bx        by  K3' 


where 

0  „   0 


U  =  \udz;V   =   \vdz  (4) 


-h  "-h 

is   the   volume  transport   and   h    =   h(x,y)    is    the   depth    of   the  water. 
§    is   assumed  small    compared   to   h. 
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Neglecting  the  effects  of  wind  stress  on  the  motion  results  in 
the  following  linearized  boundary  conditions  on  u  and  v: 

4 

<u>x=-h  "  °  »   <">2=-h  "  °  ■  <6> 

For  this  development  it  was  assumed  that  initially  a  state 

* 
of  rest  exists: 

(u)t=0  =  o  ;   (v)t=0  =  0  .  (7) 

As    in    Welander '  s    work   the   following    complex   notation  was 

introduced: 

w  =  u   +    iv 

(8) 
bn         bx  by 

Where  i  =  V^T.   Using  this  notation  eqs .  (1),  (2),  (6),  and  (7) 
can  be  expressed  as  follows: 

*    Urn   =  -  9  H  (9) 

(10) 

(ii) 
<»>z=-h  "  °  •  <12> 

Equations  (9)  through  (12)  define  the  initial/boundary  valued 
problem  that  must  be  solved  to  specify  u  and  v.   Although  (f~~)    in 
eq.  (9)  is  not  externally  specified  it  can  theoretically  be  deter- 
mined by  integrating  eq.  (3).   This  makes  it  possible  to  treat  it 
as  a  time  dependent  forcing  term  (i.e.  «j*-  (t)). 
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bw 
bt  " 

k 

*2 

b    w 

x    2 
oz 

C&.N)  =  ° 

(W)t: 

=0 

=   0 

B.       METHOD  OF    SOLUTION 

To   solve   the  problem  begin   by   multiplying    eqs .    (9)    through 

ift  *  ift 

(12)    by   e   '       and   introducing   the   relation  w      =   we         .      Then   eqs. 

(9)    through    (12)    can   be  written   as: 

I  *  2    * 

bw  .     b   w      _  ift   b?  .      . 

_     _   k  _   _   ge  g-  (t)  (13) 

QZ 

<»*>t=o  ■  °  <15» 

(»*)z=-h  -  °  •  (16> 


A  solution    to   eqs.    (13)    through    (16)    can   be    obtained  by   using 

a   time   dependent   Green's    function   approach  as    described  below. 

Assume 

w*    =   0    (t)    Z    (z)  (17) 

where  Z    (z)    satisfies    the  boundary    conditions    and  0    (t)    satisfies 

mv  mv     ' 

the    initial    conditions.      Z    (z)    is    the   eigen-f unction    determined  by 

mv  '  * 

applying  standard  separation  of  variables  to  the  homogeneous  part 
of  eq.  (13)  using  eqs.  (14)  and  (16).   It  is  found  to  be: 

Z  (z)  =  cos(pz)  (18) 

wher  e  . _   _ . 

P-  ^=111  (19) 

0    (t)    is    determined  by    substituting    eq.    (17)    into   eq.    (13) 
requiring   that   0    (t)    satisfy   eq.    (15),    and  using   the   orthogonality 
of   the   cosine  function   to    determine  the   appropriate    constants. 

Thus 

a  (tJ   =  23Lg£\     e"kP    <t-f)e-iff   6L(t.)dt,    .  (20) 

n»x  ph  J0  0nv       ' 
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Therefore  a  solution  to  eqs .  (13)  through  ( 16 )  is 

(x,y,,,t)  -^  l=i    i-'l'-IP)  ^.Vt^'l.-^^tMdf  (2!) 


w 


For    the   general    time   dependent    solution    of   eq.    (9),    dividing 
eq.    (21)    by    e  yields: 

(x.y.z.t)    =^J,    (-1  >'<">■  (P£l  fVl^UX*"*')   &f  )«■      (22) 
v     "'    '    '  h    m=l  p  J_  onv       '  v       ' 


w 


C.   DISCUSSION  OF  THE  FORMAL  SOLUTION 

It  is  apparent  at  a  glance  that  eq.  (22)  satisfies  the  boundary 

and  initial  conditions.   Substituting  eq.  (22)  into  eq.  (9)  yields 

as  a  residue: 

oo        m 

£  .  -^~) —  cos(pz)  =  -  1  .  (23) 

m=l   ph        vt-  /  \    ^/ 

It  is  easily  verified  that  the  left  hand  side  of  eq.  (23)  is 
simply  the  cosine  expansion  of  the  right  hand  side,  and  therefore 
eq.  (22),  in  the  limit,  is  an  exact  solution  to  the  initial 
problem. 

Equation  (22)  represents  the  horizontal  currents  at  all  depths 
resulting  from  the  pressure  gradient  produced  by  the  slope  of  the 
sea  surface  at  a  particular  point  and  instant  in  time.   The  integral 
term  carries  a  running  summation  of  the  currents  resulting  from  the 
slope  in  the  past  history  of  the  flow.   These  currents  are  modi- 
fied in  time  by  a  rotation  and  decay  represented  by  the  weighting 

function 

e-(kp2+if)(t-t'). 

This  weighting  function  introduces  two  interesting  time  scales 
into  the  problem.   The  first  is  the  inert ial  rotation  period,  which 

15 


varies   as   a    function    of    latitude.       For    mid-latitude    locations    it    is 
on   the  order    of   20   hours.      The  second   time   scale    is    introduced  by 
the  exponential    decay    (T    )    represented  by 

Td    =  ^  -  ^2    2    '  W 

d        kp  k(2m-l)   tt 

This    is   a   function    of   depth   and   friction,    but    for    typical    values 

2 
(k  =    .01   m  /sec,    h   =    50m)    it    is    on   the   order    of   28   hours    for    the 

first   term   in   the   series    to   decay   to   1/e    of   its    initial   value. 
This    is    the   slowest    term   to    decay. 

The  motion   represented  by    (22)    depends    on    the  surface   slope 
and   thus    new   time   scales    are    introduced   into   the  solution.      These 
depend   on    the  horizontal    dimensions    represented  and   enter    through 
the   integration   of   the   continuity    equation    (3).      Time   scales    re- 
lated to  set    up   of   the    sea   surface  and   seiche  periods    will   govern 
the   time   rate   of   change    of   the   surface  slope.      This    change 
initiates    the   transients   present    in   the   geostrophic   and  bottom 
currents,    which   show   the    inertial   rotation   and  exponential   decay 
previously    noted. 
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III.       THE   NUMERICAL    SOLUTION 

A.       FORMULATION   OF   THE  MODEL 

1 .      Establishment    of  a    Representative   Cross    Section 

In   formulating   the  numerical    scheme    it  was    decided   to 
simplify    the   model   by   neglecting    the   derivatives   along   the   y-axis. 
To  accomplish   this   without    the   loss    of  any    essential  physics,    it 
was    decided   to  allow   the   seiche   to   oscillate  across   a   narrow, 
shallow   channel    of   infinite   length.      This  prevents    sea    surface 
slopes    from   developing   along    the   length    of   the    channel    but    is 
otherwise   non-restrictive.      The  simplification    further    reduces    the 
problem  to   one   of   looking   at   a    cross-sectional    distribution   of 
variables    in   the  channel. 

It   was    further    decided   to  restrict    the  model    by    considering 
a    cross-section   of   constant   rectangular    dimensions.      The  width  and 
depth   of    the  basin    considered  were   fixed  at    5000  meters    and   30 
meters    respectively,    and   Figure    1    is   a  pictorial    depiction    of   this 
cross    section. 

These  simplifications    allow   eqs .    (3)    and    (22)    to  be   re- 
presented   in   the  model   as 

M  *  g  -  0  (35) 


,(x,z,t)    -  22  E   ,    (-1)"c°«(p»)   rV(kp2+")(t-f )  |i(f)df  (26) 

h     m=l  p  J  0xv       ' 


17 


30  m 


BASIN   CROSS-SECTION 


5000m 

z 


->x 


Z:o 


z:-h 


X=o 


MODEL  CROSS -SECT ION 


X=L 


Figure   1.      Pictorial    Depiction    of  Model    Basin 
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2 .       Establishment    of    Governing    Non-dimensional    Ratios 

To   represent    the   various    time   scales    and   frictional    co- 
efficients  possible   with   more   general    time    dependent    motion   the 
model  was    scaled  to   the  period    (T)    of  a   free    oscillating   seiche    in 
the  model    basin,    and  the   coefficient    of  vertical   friction    determined 
by   Ekraan's    "depth   of  frictional   resistance    (D)".      Equation    (26)    was 
scaled  as    follows: 

T    =  Jk  .  (27) 

fgh 

D    =   TT   jli    ;  (28) 

Let, 


t 

= 

* 
Tt 

t* 

= 

* 
Tt1 

dt' 

= 

T    dt' 

* 

* 
Lx 
x  =  —  (eleven  grid  points  across  the  channel) 

..    10A   vP* 
b£  _  o  b£ 

bx    L  "  v  * 

Ox 

*  *         *        * 

where   t    ,    tf    ,5    »    x    ,    are   non-dimensional    quantities   and  A      is 

o 

the   free   surface   elevation   at    x  =   0,    t    =  0.      Substituting    eqs .    (27) 
through    (29)    into  eq.    (26)    yields: 

w   =    2£S         (-1)    CQ5(P2)    o_  re-(-^-   +    if)(t-t')    T 

h    m=l  p  L        JQ        2tt 

be*        * 

X  -%   dt'       .  (30) 

bx 

Recalling  that  p  =  '  m~  '   ~,  and  a   =  -—-,  eq.  (30)  may  be  written 
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i„  „  rt/T      r0    ,f*     (2m-l)  ,D,    ^   . 

40g   AT00       .     1Nm         ,         f  -12tt(-:)    -^ *—  (-)    +  i 

^      o     _      (-1)    cos(pz)\  v<7 '  L      o  vh' 

e_n        /  ->_    t  \~         Ja       e 


L  m=l      (2iq-1)tt        ^0 


2       _    2 


(t-t»)* 


* 

x  ^  dt'    .  (31) 

bx 

It   was    determined   that    the   two   non    dimensional    ratios, 
(f/CT)    and    (D/h),    appearing    in   eq.    (30)    govern    the    time/spacial 
scale  and  frictional   dependence   of    the   motion.       For    discussion   and 
computation    the    first   was    designated    (PCF)    and   the   second    (PCD)  . 

The   first    of   these   ratios    represents    the   relationship 
between    the   earth's    inertial    frequency    and   the   frequency    of   the 
basin    in    the   following   manner . 

f   =   2fisin   0 

where  £2  is    the   frequency    of   the   earth's    rotation   and   0   is    the 

latitude. 

a    =   2tt  F    , 
m 

where  F      is    the  resonance  frequency    of   the  basin.      Since  no  external 
m 

forcing    terms   were  used    in   the   solution,    F      is    also   the   approximate 

m  rr 

frequency    of   the  motion.      Thus: 

f   _  Q  sin0 

OF 
m 

The  model  was  scaled  for  (0  £  PCF  ^  1).   A  value  of 
(PCF  =  0)  represents  motion  of  a  time  scale  which  occurs  too 
quickly  to  be  affected  by  Coriolis  accelerations,  while  a  value 
of  (PCF  =  1)  represents  motion  whose  time  scale  is  equal  to  that  of 
inertial  rotation  and  subsequently  greatly  influenced. 


20 


The  second   term,    PCD,    is    the  ratio  of   the    depth  of 
frictional    influence    (Sverdrup,    Johnson,    and   Fleming,    1942)    to   the 
actual   water    depth.      The  coefficient    of   vertical    friction  was    de- 
termined  in   the  model    from   eq.    (28)    as 

k  =   (PCDX?h)2   Xf    .  (3la) 

2n 

The  model  was  scaled  for  (0  ^  PCD  ^2).   A  value  of  (PCD  =  0) 
represents  a  very  deep  basin  where  frictional  effects  on  the  motion 
are  minimal,  while  (PCD  =  2)  represents  a  basin  where  friction  domi- 
nates the  motion. 

If  the  time  scale  of  the  motion  being  represented  by  the 
model  is  controlled  strictly  by  the  resonance  of  the  basin,  then 
the  spacial  scale  of  the  basin  is  determined  by  a  choice  of  PCD 
and  PCF.   PCD  will  govern  the  approximate  depth  of  the  basin  while 
PCF  determines  the  basin  width  through  eq.  (27).   Certain  con- 
clusions about  motion  caused  by  external  forces  (i.e.  tides, 
atmospheric  pressure)  producing  time  scales  independent  of  the 
basin  size  were  also  inferred  from  the  model  and  are  discussed  in 
the  conclusions. 

3«   Adaptation  of  the  Formal  Solution  to  the  Model 

To  establish  two  equations  that  could  be  solved  simulta- 
neously for  sea  surface  elevation  and  transport,  eq.  (26)  was 
integrated  in  the  vertical  using  eqs .  (4)  and  (8).   Thus: 

p  ° 

and  W  =  U  +  iV.   Equation  (32)  represents  the  total  volume  trans- 
port associated  with  the  flow  at  a  particular  point  in  the  channel. 

21 


The  real   part    of   eq.    (32)    represents    the   cross-channel    transport, 
while   the    imaginary   part    the  along-channel    transport. 

The   numerical    scheme  was    established   to  solve   eqs .    (25) 
and   (32)    for    §   and  W.      The   velocity  profile,    bottom   stress,    and 
energy   balance  associated  with   the   flow   were   also  recoverable   from 
the   scheme,    and  the  details    for    this    recovery   will   be   discussed 
later . 

4.  Initial    Conditions 

The   model   assumes   an    initial    state   of  rest    in    the   motion 

and  an    initial   set-up    of   the  free   surface.      This   was    represented 

as 

W(x,0)    =  0,  (33) 

and 

?(x,0)    =  Aocos  &Zl   }  (34) 

where  W   =  W(x,t)    and  §    =   ?(x,t). 

5.  Boundary   Conditions 

Since  the  governing  equations  were  solved  without  reference 
to  horizontal  boundaries,  these  conditions  had  to  be  built  into  the 
model.   The  model  therefore  assumes  no  flow  across  the  horizontal 
boundaries,  or 

W(0,t)  =  W(L,t)  =  0.  (35) 

These  conditions  were  satisfied  by  requiring  the  surface  slope  to 
remain  zero  at  these  boundaries  at  all  times.  Therefore  eq.  (32) 
is  forced  to  meet  the  conditions  of  eq.  (35). 

6 .  The  Step-by-step  Scheme 

Both  Welander  and  Gait  suggested  that  a  step-by-step 
numerical  scheme  could  be  used  to  solve  eqs.  (25)  and  (32).   As  a 
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first    solution    the   staggered   central   difference   scheme   suggested 
by   Gait   was    used   and    is    summarized  below. 

t    =    -    At        Given    initial    conditions    on   W. 

t    =    0  Given    initial    conditions    on   §. 

t    =    At  Calculate  W   from  eq.    (32)    using  r*  at   t   =  0, 

t    =   2At  Calculate  £    using  an   explicit    finite-difference 

scheme   on    eq.     (25)    and   U   at    t    =   At.       Calculate   a    new 
surface  slope. 
t    =   3At  .    .     .    etc. 

The   term   At    is    the   time  step    of    the  explicit   scheme. 

This    scheme  was    developed   neglecting   Coriolis    (f    =   0)    in 
eq.    (32).      The   scheme  proved  to  be   stable  and  accurate   subject   to 
the   following    stability   criterion: 

YIh^<   1    •  (36) 

The   expression    Ax    is    the   spacial   step   across    the   channel.       It   was 
felt    however    that   a  more  precise   solution,    not   subject    to   such  a 
restricted  stability   criterion,    could  be   formulated  using   a    totally 
implicit    finite-difference   scheme.      The   development    of   this    method 
is   now  presented. 

7.      The    Implicit    Finite-difference    Scheme 

The    implicit   scheme   used   is    a    modification    of  a   method 
described  by   Richtmyer    and  Morton    (1967). 

Start   by  approximating   the   time  rate   of   change   of   the 
transport   as    follows: 


bt  At 


(37) 
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From   eq.    (32) 


bw 
bt 


P      L^O 


=  I.  ("2fl  s       JL 

At    /     h    m=l      : 

L  P 

^  e-(kp2+if)(t-tM  |L(t.)dt, 


t+At  e-(kp2+if)(t+At-f )   6S(t.)dt 


(38) 


By   expanding   the   first    integral   as   \      +  \  and  regrouping    terms 

J0       Jt 
eq.    (38)    may   be  written   as 


»S  =  JL  f-22  e       1   C*  e-(Kp2*«)(t-f)(e-(kp2-if)At.     M.(t,)dt,\ 

ot         At   I      h    m=l      2   J.  v  '    0xv       '         J 


At   (     h     m=l    2  ^t 


At   e-(kp2+if)(t+At-t')    b£ 


bxx       ' 


(39) 


Eqs .    (25)    and    (32)    can    then   be   represented    in    the  following   form: 


^bt;       At  2 


bu\         /bu 
bx;t       l^bx 


(40) 


t  +  At  J 


where 


i$3L  =  #S_%  +  /C2_ 

vbt;         At         vAt     ;        v2At' 
t      2 


2     . 


fill  +/ss 

bxjt        ^6x, 


2    . 


(41) 


a».-afi.  1  f  .-civ  +if)(t-t»>  -(kp+if)At     s(t.)dt.  (42) 

h    m-1      2  J,.  v  '    0xv       ' 

P        ° 


00  x        pt+At 

=      m=l      ~2      \ 
P  t 


C2*    =    "42s    ,      i     V" """   e-(kp"+if)(t+At-t')dt, 


24 


Equations  (40)  and  ( 4l )  can  be  represented  in  Matrix  form  as 
follows : 


<t?'t+At=A  +  f^  [(6)t  +   <eW]  (*3) 


where, 


0   = 


B   = 


0      0      C2/At 
0      0      D2/At 
■10  0 


and, 

CI    =   Real    {ci    }     ;       C2    =   Real    { 02    } 
Dl    =    Imag    {ci    }     ;       D2    =    Imag    { C2    } 

(Real   {    }    and    Imag   {    }    represent    the  real   and    imaginary  part   of 
the  argument    in    braces). 

The   following    implicit    finite-difference   equation  was    used 
to  represent    eq.    (43): 


0n+1-   0n 


D_   _ 


At 


=   A   + 


B 

4ax 


9n  +  *    -    9n+^    +   9n    ,     -    9n    , 
j+1  j-1  3+l  j-1 


(44) 


where  9 . 
3 


QnAt 
JAx 
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and  n  =  (1,2,3, •••,  maximum  number  of  time  steps) 

j    =  (1,2,3,...,  J) 

are   the  grid  points  in    the  time  and  x-directions   respectively. 

Letting 

a.    =  -, ,    and 

4Ax 


D.    =  tt9n    ,     +   9n 


-   <*9.    _    +    AtA, 


eq.    (44)    becomes 


(45) 


The    left   hand   side   of   eq.    (45)    contains    the  unknown   values    of 
(9)    to  be   calculated,    whereas    the   right   hand  side    (D.)    are  all 
known   values. 

A  recursion   relation    of   the   following    form  was    assumed. 


9  .         =   E.   9  +    F  . 

D  J      3+1  3 


(46) 


where 


E.   = 


eJ        eJ 
ell      e12 


'21 
J 


'22 
J 


13 
'23 


31        32        33 


and 


F.    = 
3 
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If  eq.    (46)    holds    for   all   9.,    then    the  model's    left   hand  boundary 
conditions   require   that 


0 

0 

0 

E      = 
o 

0 

0 

0 

0 

0 

1 

and 


F      = 
o 


(47) 


,n+l 


Substituting   the   expression   for    9.        derived  from   eq.    (46)    into 
eq.    (45)    yields: 

D.-a    F 


9 


n+1 


a 


i+a  e.  , 


±1 


i+a  e 


J'1  .. 


(48) 


Equating  the  right  hand  side  of  eq.  (48)  to  eq.  (46)  gives  the 
following  recursion  relationship: 


E.  =  M 1 

3    I1**  Ej-i. 


D.-Q!    F.      ' 
-2 111 


Fj 


j-l    J 


,     j  *  i; 


>     j  a  i; 


(49) 
(50) 


and   I    is    the    (    3X3    )    identity   matrix. 

Using   eqs .    (47),    (49),    and   (50),    values    of   E.   and  F.   were 
calculated   inductively    in    order    of    increasing   j(j    =    1,2,3,...,    J-l). 

Since   the  value   of  9 .    ,    is    given   for    j    =   J-l   by    the  model's    right 

j+1  J 

hand  boundary   condition    (U      =   V      =   0   and   %,  -   %       ),    values    of 

9.        were   then   calculated   inductively    from   eq.    (46)    in    order    of 
decreasing    j(j    =   J-l,    J-2,     ...     ,    1).      This    completes    the   calculation 

y.  i        1 

of  9 .    (j  =  1,2,3,  ...  ,  J-l).   The  details  of  the  mathematics  for 
representing  E.,  F.,  D.,  and  eq.  (46)  are  presented  in  appendix  A. 
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8.      Numerical   Treatment    of  CI      and  C2 

The   numerical    scheme   for    representing   CI      and  C2      is 
extremely    lengthy   and    is    merely    outlined   in    this    section.      For    a 
more   detailed  presentation    the   reader    is    referred   to  appendix   B. 

The    integration   of  CI      was  performed  by    dividing   the   total 
integration   time    into   small    time    intervals    (At)    and  assuming    that 
the   surface   slope  and  effects    of   inertial    rotation    could  be   re- 
presented adequately   by  a    constant   throughout  ^:he    interval.      These 
values    were  calculated  at   the   raid-point    of   each   At    interval  and 
used   in    the   scheme  as    representative   values    throughout    the  entire 
sub-interval.      This   permitted   the    integral    to  be   evaluated  analyt- 
ically   over    each    interval   with   a   maximum   error    on   the   order    of 
these  approximations.      Using   these  assumptions    CI      was    numerically 
represented  as 


cl*i  =  -&  rv-L 


j_(1.e-  kP2At) 


M=l  (^p      Ikp  J  N=l  L 


i-i  r  ^B   ..  i  .2 

S 


,^^-h,    -    kp    At    -if At   .. 


2 


X   (e"if(I"^"N)At)(e'kp    (I-1-N)At 


(51) 


*1 

for    I    =   2,    3,    ...    ,   and  CI        =0.      The   value   of   the    integer    I    ranges 

from   1    to  MAXT    (maximum  problem   time)    and   IAt    represents    the   elapsed 
time   from    initial    conditions.      L   takes    on    integer    values    from   1   to 
MAXL    (maximum   number    of    cross-channel   grid  points)    and  LAx    is    the 
distance   from   the   left   boundary.      The    integer    M   ranges    from   1   to 
MAXM,    and   determines    the   maximum   number    of  modes    in    the   series    that 
are    included   in   the  approximation   of   CI    .      It   was    determined   from 
eq.    (23)    that   MAXM   =   25   would  represent    the    solution  with   a   maximum 
error    of    three  percent,    and   this    value   was   used  throughout   all 
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numerical  evaluations.   The  integer  N  takes  on  the  values  from  1 
to  (1-1)  and  determines,  through  the  exponential  terms,  how  much 
rotation  and  decay  has  taken  place  in  the  recent  time  history  of 
each  transient  considered.   Since  the  terms  of  the  series  decay 
more  rapidly  with  increasing  friction  (especially  the  higher  modes), 
it  was  not  necessary  to  carry  their  integration  as  far  back  as 
N  =  1.   It  was  decided  to  include  only  the  terms  of  significant 

magnitude,  and  this  was  accomplished  by  requiring  that  when  the 

2 
value  of  kp  (I-l-N)At  in  eq.  (51)  was  greater  than  10,  the  cor- 
responding terms  were  discarded. 

The  integration  of  C2   involves  only  one  time  interval 
(t  to  t+At) .   Again  the  amount  of  inertial  rotation  was  evaluated 
at  the  mid-point  and  assumed  constant  throughout  the  interval,  and 
02      was  evaluated  as  the  following  constant: 


C2 


_   MAXM    .  ,.,!  A„_x 
•  =  - 2g_  s     e-if(?§At) 

h  M=l 


2-  (i-e-kP2At) 


kp' 


(52) 


CI   and  C2   both  contain  the  term 


— —  (1-e  p    )   ,  and  for  k  =  0  this  term  can 
.  kp  J 


kp 
not  be  numerically  evaluated  in  this  form.   However  this  difficulty 

2 

.         , .     -kp  At  .     _ 
was  overcome  by  expanding  e  c  in  a  Taylor  series,  i.e. 


(i  -  kp2At  *  i^n\  us^Ati3) 


(53) 


and   substituting    this    series    for    the   term    in   eqs .    (51)    and    (52) 

2  -U 

whenever    the   value   of   kp    At   was    less    than    10      . 

The   method   for    separating    the  real    and    imaginary   parts    of 

*  * 

CI      and  C2      involves    using   the   relation 

e  =   cos(ft)    -    i   sin(ft) 
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in   eqs .    (51)    and    (52)    and   gathering   terras.      The   details   are 
covered    in   appendix    B. 

9.  Velocity  Calculations 

The  velocity  profiles  across  the  channel  were  calculated 
by  the  model  at  discrete  time  intervals  determined  by  an  input 
parameter.   Noting  the  simularity  between  the  integral  of  eq.  (26) 
and  CI  ,  it  was  numerically  possible  to  compute  and  store  the  terms 

of  these  integrals  needed  for  velocity,  during  the  calculation  of 

* 
CI  .   For  velocity  profile  computation,  these  stored  values  were 

recalled  and  applied  in  eq.  (26).   The  details  for  numerically  re- 
presenting eq.  (26)  is  also  presented  in  appendix  B. 

10.  Bottom   Stress    Calculations 

The  bottom  stress  (T)  associated  with  the  motion  was 
calculated  using  the  following  relation: 

'  = k  (If )  ...* 

or  fr om  eq.  (26) , 


r=Ms"   f  e-(kP2+if)(t-t-)  jji(t.)dt.  . 

h   m=l  Jn  bxv  ' 


(54) 


Since  the  integral  term  of  eq.  (54)  is  identical  to  that 
of  eq.  (26),  bottom  stress  is  easily  calculated  whenever  the 
velocity  profile  is  recovered.   Appendix  B  also  contains  the 
details  for  this  calculation. 

11.   Energy  Calculations  , 

The  total  energy  balance  associated  with  the  seiche  is  a 
sum  of  its  potential  and  kinetic  energy.   This  balance  was 
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calculated  as  follows: 


PE  =  h  p   g  \      %2   dx  (55) 

J0 

KE  =  h   P  \    \   w   dx  dz  ( 56 ) 


PE  and  KE  are   the  potential   and  kinetic    energies   and    (p)    is    the 

is   the  density    of   the  water.      Equations    (55)    and    (56)    were   evaluated 

using   Simpson's    1/3   rule    (McCalla,    1967)    whenever    velocity  was 

calculated. 

The   energy    calculations   were   monitored   during  model 
testing   as   a   calculational    check  since   small   numerical   errors 
were   quickly   apparent    in    increasing    energy    levels. 
12.       Scaling   of  Model  Output 

Since   the   governing    equations    for    the  model   were   scaled 
to  produce   the   ratios   PCF  and   PCD,    the   output    of    the   model   must 
be  appropriately    re-scaled   to   obtain   real    values    for    the   variables. 
The  series    of  graphs    were  produced   on   a    CALCOMP   plotter    with    cor- 
responding  values    from   the  model.      The   various   axes    were  auto- 
matically  scaled  by    the  computer    to   fit    the   data   being   represented. 
These    axes     were    then   re-scaled  using   a    non-dimensional   scaling 
factor    derived  from   the   continuity    equation   to  allow   real   values 
to  be   easily   obtained.      The   exception    is    the  figures    depicting   the 
energy   balance,    which  are   not    scaled  and   should   only    be   used  to 
indicate   relative    changes    in    energy    levels. 
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The  scales    for    the   various   axes   are  as    follows: 

AXIS  SCALE 

Surface   Elevation  §/A 

o 


Time 


Velocity 


t    (gh)2 


2L 


2w    (h/g) 2 

C..A 
1    o 


2W 
Transport  1 — 

C2Ao(gh)^ 

2T(h)3/2 
Bottom  Stress  '   i 

C3Aokg^ 

The  values    of  g ,    h,    k,    A    ,    L  are   determined   from   the 
actual   basin   being    represented   and  C.  ,    C    ,    C     are   scaling   con- 

J-  £•  J 

stants  dependent  on  the  parameters  of  the  prototype  basin  and  the 
scaling  of  the  graphs  by  the  computer.   These  values  are  calculated 
for  each  figure  and  inserted  in  the  scale.   The  variables,  t,  w,  W, 
§,  and  T  are  the  real  values  associated  with  the  basin  being  re- 
presented and  may  easily  be  obtained  from  the  graphs  by  entering 
appropriate  arguments  for  g,  k,  h,  A  ,  and  L  in  the  corresponding 
scale. 

B.   TESTING  THE  MODEL 

1.   Establishing  a  Reference  Solution 

One  of  the  major  difficulties  in  numerical  analysis  is 
verifying  the  accuracy  and  stability  of  the  scheme.   Since  many 
numerical  solutions  deal  with  problems  that  have  no  analytical 
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solutions  it  is  common  practice  to  test  the  accuracy  of  these 
models  against  empirical  data,  or  analytical  solutions  involving 
degenerate  cases.   As  previously  noted,  data  pertaining  to  time 
dependent  flow  is  extremely  scarce,  and  it  was  decided  to  test  the 
accuracy  of  the  model  against  a  known  analytical  solution  for  the 
degenerate  case  where  PCD  =  0  (i.e.  no  friction).   Analytical 
solutions  to  this  case  are  readily  available  in  the  literature  and 
a  modification  of  one  proposed  by  Prodman  (1952)  was  used  in  the 
verification  of  the  accuracy  of  the  numerical  scheme. 
2 .   Accuracy  and  Stability  of  the  Model 

The  accuracy  of  the  model  was  determined  by  comparing  the 
numerical  representation  of  the  free  surface  elevation  at  (x  =  0), 
the  raid-channel  transport,  and  the  total  energy  balance,  with  the 
analytical  values.   Figures  2  through  5  highlite  this  comparison 
and  indicate  that  the  accuracy  of  the  model  was  a  function  of  time 
step  (At)  and  the  Coriolis  parameter  (f).   It  was  noted  that  as  the 
time  step  was  reduced  the  numerical  solution  appeared  to  converge 
to  the  analytical  solution.   The  phase  shift  in  surface  elevation 
over  five  periods  was  reduced  from  45  degrees  for  (At  =  T/20)  to 
24  degrees  for  (At  =  T/40)  and  appeared  to  be  independent  of  the 
value  of  Coriolis.   This  error  might  possibly  be  introduced  into 
the  system  by  approximating  the  sea  surface  slope  as  a  constant 
value  in  the  time  history  integration  scheme.   As  (At)  increases 
this  linear  approximation  is  extended  over  a  greater  time  span 
and  the  resulting  errors  amplified.   As  the  value  of  the  Coriolis 
parameter  was  increased,  the  error  in  representing  the  amplitude 
of  the  free  surface  and  total  energy  balance  also  increased,  as 
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depicted  in  Figures  3  and  5-   The  amplitude  deviation  was  a  maximum 
of  one  percentage  for  (At  =  T/20)  and  the  energy  rise  was  8.5  and 
3.7  percentage  for  (At  =  T/20)  and  At  =  T/40)  respectively.   Errors 
of  this  nature  did  not  appear  for  solutions  which  neglected  Coriolis, 
and  possibly  were  introduced  by  approximating  the  amount  of  rotation 
[i.e.  e    *     ' ]  as  a  constant  in  the  integration  scheme.   It  is 
not  exactly  clear  how  this  error  increases  the  energy  in  the  model 
but  it  was  observed  that  it  appeared  linear  and  very  systematic. 

To  determine  whether  this  energy  rise  constituted  a 
numerical  instability,  a  standard  Fourier  series  method  (Smith, 
1965)  of  stability  analysis  was  attempted  on  eq.  (44).   However, 
due  to  the  complicated  nature  of  the  amplification  matrix,  the 
eigenvalues  could  not  be  readily  obtained  and  the  analysis  was 
abandoned.   As  an  alternate  test,  the  model  was  run  with  (At  =  T/5) 
for  (t  =  5T) .   This  corresponds  to   fgh  —  >  3   ,  and  resulted  in  a 

l      AX       _t 

further  linear  rise  in  the  energy  level.   These  errors  were  there- 
fore judged  to  constitute  a  problem  in  numerical  accuracy  and  not 
stability,  for  the  range  in  which  the  test  cases  were  run. 

It  was  further  found  that  the  magnitudes  of  these  errors 
could  be  predicted  and  controlled  by  varying  the  time  step,  and  the 
model  was  run  maintaining  a  maximum  allowable  error  of  one  percent 
in  the  elevation  and  transport  amplitude,  and  four  percent  rise  in 
the  total  energy  balance.   This  was  accomplished  by  using  a  time 
step  of  (At  =  T/20)  for  (PCF  <  .5),  and  (At  =  T/40)  for  (PCF  ^  .5). 
Since  the  errors  in  phase  shift  were  strictly  a  function  of  time 
step,  they  were  completely  predictable  and  accounted  for  during 
data  analysis. 
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It  should  be  noted  that  for  cases  where  friction  is 
neglected  (k  =  0),  the  transients  do  not  decay  with  time  and  must 
be  carried  throughout  the  integration  scheme.   Additionally,  the 
bottom  boundary  condition  degenerates  to  a  free  slip  situation  and 
the  approximation  of  a  constant  velocity  profile  with  a  truncated 
cosine  series  becomes  more  inaccurate.   The  number  of  calculations 
increase  rapidly  as  friction  is  decreased,  which  may  introduce 
increasing  errors  in  round-off. 

Based  on  the  comparison  with  the  analytical  solution  it 
was  concluded  that  the  accuracy  of  the  model  could  be  sufficiently 
controlled  such  that  the  details  of  time  dependent  motion,  as 
described  by  the  governing  integral  and  differential  equations, 
could  be  adequately  represented  and  examined. 

C.   MODEL  LIMITATIONS 

The  most  serious  limitation  imposed  on  the  model  were  the 
basic  assumptions  of  Ekman  dynamics  used  in  the  formal  solution 
of  the  governing  equations.   These  assumptions  were  considered 
necessary  to  simplify  the  mathematical  treatment  of  these  equations 
to  the  extent  that  a  numerical  solution  was  feasible.   The  additioanl 
simplifications  imposed  during  the  formulation  of  the  model  further 
limit  its  capabilities.   It  was  felt  however  that  valuable  experi- 
ence in  the  techniques  of  numerically  representing  time  dependent 
motion  would  be  gained  during  the  formulation  of  this  simplified 
model,  and  insight  gained  from  the  numerical  product.   With  this 
information  available  it  might  be  possible  to  further  refine  the 
model  and  remove  many  of  the  more  important  restrictions.   The 
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limitations    inherent    in    the  present   model   are   listed    in   this 
section   and   the   more    important    ones    will   be   discussed    in   section 
(V)    with   recommendations    for    improving   the   scheme. 

The  basic  assumptions    made    in  arriving   at    eq.    (22)    have   the 
following    limiting    effects    of   the  model's    capabilities:       1)    the 
non-linear    terms    in   the  motion   have   been    ignored;       2)    a    variation 
in   the   coefficient    of   vertical    friction    in    the   layer    was    not 
permitted;       3)    the   assumption    of    constant    density    neglects    the 
effects    of   stratification   on    the  pressure   gradients    in   the   layer, 
and   the   subsequent   modification    of   the   vertical   momentum   transfer; 
and      4)    horizontal    friction   was    considered    insignificant   and 
neglected. 

Simplification    of   the   model    resulted    in    the  additional 
limitation:       5)    the  assumption   of  a    constant   basin   neglects    the 
important    effects    of  bottom   topography    on    the   flow.      This    is    es- 
pecially   important    for    considering   motion    in    shallow   near    shore 
areas   where   these   effects   are   large. 
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IV.   PRESENTATION  OF  RESULTS 

A.   TEST  CASES  AND  GEOPHYSICAL  SIMULATIONS 

1.  Description  of  Test  Cases 

In  order  to  gain  maximum  insight  into  the  physical  processes 
occurring  in  the  motion,  as  well  as  to  study  the  effects  of  inertial 
rotation  and  friction  on  the  variables,  the  model  was  run  for  25 
test  cases  in  the  ranges  of  (0  £  PCF  ^  1)  and  (0  £  PCD  ^  2). 
Table  I  summarized  the  various  controlling  parameters  of  each  run 
and  includes  the  Central  Process  Unit  (CPU)  time  used  on  the  IBM 
36O/67  system  for  each  program. 

2.  Description  of  Geophysical  Simulations 

As  interesting  dynamic  interactions  developed  in  the  test 
cases,  it  was  decided  to  simulate  actual  geophysical  embayments 
to  determine  to  what  extent  these  phenomena  were  present  and 
significant.   The  actual  length,  depth,  and  latitude  of  the  basins 
were  measured  and  used  to  calculate  (f )  ,  (0-),  and  (PCF).   An  Ekman 
depth  of  20  meters  was  assumed  in  the  calculation  of  (PCD) .   Runs 
26-3O  represent  these  simulations  and  the  various  imput  parameters 
are  summarized  in  Table  II. 

B.   DISCUSSION  OF  TEST  RUN  RESULTS 

1.   Frictional  Effects  on  the  Flow 

To  describe  the  conditions  where  friction  totally  dominates 
the  motion,  the  model  was  run  with  varying  amounts  of  friction  and 
(PCF  =0).   It  must  be  noted  that  by  neglecting  the  effects  of 
inertial  rotation,  the  Ekman  depth  is  no  longer  specified  as 
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Table   I.      Parameters    of  Test   Cases    (runs    1-25) 


RUN   NO. 

PCF 

PCD 

CPU(MIN-SEC) 

1 

.25 

.5 

6-2 

2 

.25 

1.0 

3-51 

3 

.25 

1.5 

3-8 

4 

.25 

2.0 

2-49 

5 

0 

-5 

5-H 

6 

0 

1.0 

3-20 

7 

0 

1.5 

2-4? 

8 

0 

2.0 

2-24 

9 

0 

0 

13-33 

10 

.25 

0 

14-58 

11 

.5 

2.0 

6-25 

19 

.75 

2.0 

6-0 

13 

1.0 

2.0 

5-48 

14 

.5 

1.5 

7-25 

15 

.75 

1.5 

6-42 

16 

1.0 

1.5 

6-25 

17 

.5 

1.0 

9-40 

18 

.75 

1.0 

8-21 

1Q 

1.0 

1.0 

7-40 

20 

•  5 

.5 

16-14 

21 

.75 

.5 

13-52 

22 

1.0 

.5 

12-29 

23 

.5 

0 

62-48 

24 

.75 

0 

62-38 

25 

1.0 

0 

62-51 

NOTES:    For    the  model   basin 
10        T   =    583.2   Sec 


2) 


CT    =    .108   X   10_1    Sec"1 
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evidenced  by  eq.  (28).  Therefore,  the  frictional  coefficients 
generated  by  runs  1  through  4  were  artificially  imposed  on  the 
model  with  (PCF  =  0),  and  resulted  in  runs  5  through  8. 

The  results  for  these  cases  are  depicted  in  Figures  6 
through  8  and  are  simply  a  seiche  damped  by  friction.   As  the 
friction  coefficient  increased  the  damping  of  the  motion  increased 
with  a  resulting  increase  in  the  period  of  oscillation.   For  the 
maximum  case  of  (PCD  =  2),  the  seiche  was  completely  damped  after 
4^2  periods.   These  results  were  completely  predictable  and  are 
presented  to  show  that  the  model's  representation  of  motion  which 
is  dominated  by  friction  conforms  to  previous  observations. 
(Pierson  and  Neuman,  1966) 

An  interesting  coupling  between  bottom  stress  and  transport 
developed  during  these  cases  and  a  typical  example  is  depicted  in 
Figure  9.   In  all  cases  examined  a  phase  lag  quickly  developed  be- 
tween the  transport  and  the  bottom  stress.   For  the  case  showed, 
the  value  of  this  lag  averaged  45  degrees  through  the  periods  of 
oscillation.   This  phase  lag  can  be  explained  by  recalling  that 
the  transport  is  an  integration  of  the  motion  over  the  entire  water 
column,  whereas  the  bottom  stress  is  solely  dependent  on  the 
velocity  gradient  at  the  bottom.   Since  friction  has  its  greatest 
effect  at  the  bottom,  changes  in  magnitude  and  direction  of  the 
flow  nearly  always  initially  occur  next  to  this  boundary  and  prop- 
agate upwards  in  the  column.   The  result  is  that  changes  in  trans- 
port always  lag  changes  in  the  bottom  velocity  gradient  and 
consequently  the  bottom  stress.   This  phase  lag  has  obvious  effects 
on  the  accuracy  of  commonly  used  numerical  schemes  where  bottom 
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stress  is  parameterized  in  terms  of  the  transport.   This  technique 
places  strong  restrictions  on  possible  interactions  between  the 
flow  and  the  bottom,  and  although  it  may  be  sufficient  for  approx- 
imating steady  state  problems,  it  is  unlikely  that  it  will  represent 
transient  motion  very  well. 

Figure  10  represents  typical  velocity  profiles  for  motion 
which  is  totally  dominated  by  the  effects  of  vertical  friction. 
The  increased  period  of  oscillation  and  the  mechanism  causing  the 
phase  lag  between  the  transport  and  the  bottom  stress  is  clearly 
evident.   At  times  (t  =  T/2)  and  (t  =  3T/4)  the  profiles  are  in 
transition  from  being  totally  positive  to  totally  negative,  and  the 
transport  and  bottom  stress  are  180  degrees  out  of  phase. 

These  velocity  profiles  are  typical  of  those  found  in  other 
boundary  layer  flow  and  are  similar  to  the  wind  profiles  observed 
at  the  air  sea  interface  over  smooth  water.   They  indicate  that  the 
stability  of  the  flow  has  a  strong  frictional  dependence.   At 
(t  =  T/4)  and  (t  =  T)  the  profiles  are  maximum  and,  except  initially 
in  the  upper  layers,  friction  exerts  a  strong  influence  throughout 
the  column.   The  result  is  a  smooth  regular  flow  and  a  well  de- 
veloped frictional  boundary  layer  with  the  region  of  greatest  shear 
at  the  bottom.   During  the  transition  phases  (t  =  T/2  to  t  =  3T/4) 
the  profile  changes  direction  and  the  water  column  appears  to  be 
less  stable.   The  velocities  in  the  column  are  minimum  just  prior 
to  and  after  reversal  and  thus  the  influence  of  friction  in  the 
column  is  reduced.   This  region  of  reduced  frictional  influence 
moves  up  the  column  at  the  point  where  the  flow  changes  direction. 
As  the  flow  completes  its  reversal  and  the  velocities  in  the  column 
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Figure    10.    Effects    of   Friction    on   Mid-channel   Velocity 
(PCF   =   0,    PCD  =   2) . 
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increase,  the  effect  of  strong  frictional  influences  is  quickly 
re-established  throughout  the  layer,  the  region  of  maximum  shear 
shifts  to  the  bottom,  and  smooth  flow  once  again  prevails  at  time 
(t  =  T). 

2.   Coriolis  Effects  on  the  Flow 

To  gain  insight  into  how  inertial  rotation  modifies  the 
transient  motion,  runs  10,  and  23  through  25  were  made  neglecting 
friction  (PCD  =  0) . 

It  is  important  to  realize  that  the  relationship  between 
the  transport  and  the  free  surface  elevation  for  a  non  rotational, 
frictionless  seiche,  as  seen  in  Figures  6  and  7,  is  similar  to  the 
velocity/elevation  relationship  of  the  frictionless  pendulum.   At 
(t  =  0),  (t  =  T/2) ,  and  (t  =  T)  ,  the  energy  of  the  seiche  is  en- 
tirely potential,  whereas  at  (t  =  T/4)  and  (t  =  3T/4)  it  is 
entirely  kinetic.   At  intermediate  times  the  energy  balance  is  a 
combination  of  potential  and  kinetic  and  is  totally  conserved.   The 
free  surface  elevation  and  transport  periodically  oscillate  between 
maximum  positive  and  negative  values  of  respectively  equal  magni- 
tudes, always  maintaining  their  constant  phase  shift.   Figure  2 
represents  this  oscillation  for  the  free  surface. 

As  the  effects  of  Coriolis  acceleration  were  introduced 
into  the  model  these  conservative  oscillations  were  greatly  modified 
as  can  be  seen  from  a  comparison  of  Figures  2,  3,  11,  and  12.   The 
effect  of  a  moving  earth  is  that  the  total  potential  energy  of  the 
surface  elevation  will  not  be  available  for  transfer  to  kinetic 
energy  since  the  water  particles  are  forced  to  turn  due  to  the 
Coriolis  Acceleration.   Although  the  energy  balance  still  remains 
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conservative,    the  relationship   between   the   free   surface   elevation 
and  the    flow   changes    significantly. 

To  understand  how  rotation   acts    to  modify    the    flow,    con- 
sider   first    the  physical   processes    which    drive   the   system.      These 
are   discussed   in    terms    of   the   frictionless    irrotational   seiche. 
The   initial    conditions   for    the  model   were   identical    for   all   runs; 
that    of  a  positive   set-up    of   the   free   surface  at    the    left   boundary 
(x   -   0)    and  no  motion    in   the   basin    (a    condition    of   static   equili- 
brium).     As    the   system   is    released  the  pressure   gradients    developed 
by   this    set-up    start    the   fluid  moving    in    the  positive  x-direction 
resulting    in   a    net    divergence    in    the    left    side    of    the  channel    and 
convergence   in   the  right.      The   free   surface  responds    to  this    flow 
by   slumping   at    the    left   and   rising   at    the   right   with   equal   magni- 
tudes.     At    time    (t    =  T/4)    the   surface    is    level    (i.e.    no  gradients) 
and   the   kinetic   energy   associated  with    the   flow    is    equal   to   the 
potential    energy    of   the    initial    free   surface   set-up.      Therefore   the 
surface   continues    to   fall   at    the   left   and  rise  at    the   right   until 
this   kinetic   energy   has    been    entirely    converted    into  potential 
energy   and  a    state    of   momentary    equilibrium    is    once   again   reached. 
The   sea    surface    is    now  set    up   at    the   right   boundary    equal    in   magni- 
tude  to    the   initial    conditions,    the  pressure   gradients    equal    but 
opposite   to  the    initial   values,    and  the    flow   reverses    direction 
and   the   system  returns    to    initial    conditions,    etc.       It    is    important 
to  realize   that    only    the   divergence   of   the  flow    in    the  cross-channel 
direction    can    influence   the   system  because   no  gradients   are  allowed 
to   develop   along    the    length    of   the  channel. 
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With   these  processes    in   mind   the   observed   effects    of   ro- 
tation   on    the   system   can   adequately    be   explained.      With   the    initial 
set-up  previously    described,    the   flow   will   start    in   the   same  manner. 
As    soon   as    the  motion   begins    however,    the  Coriolis   acceleration  will 
start    to   turn   the   flow    in    the  negative   y-direction    (left   rotation), 
and  a    component    of   the   flow  will   start    to  develop   along   the  axis    of 
the   channel.      Therefore   some   of   the  potential    energy    is    converted 
into   rotational   kinetic    energy.      The   amount    of    rotational    kinetic 
energy    that    the   system   receives    is    a    function    of    the   velocity    of 
the   flow,    the   value   of  Coriolis,    and    the   time  scale   of    the   motion. 
The   slumping    of   the  surface   elevation   at    the    left   and  rising   at    the 
right   will    continue  as    long   as    there    is    divergence    in   the  positive 
x-direction;    however    since   a    certain   amount    of   energy    is    rotational, 
the   set-up   at    the  right    boundary    can   not   achieve   the  magnitude   of 
initial    conditions.      When    the   cross-channel    flow   vanishes    the   en- 
tire  motion    is    in    the  along    channel    direction,    with   subsequent 
rotation    turning    the   flow   back    in    the    direction    from  which    it 
initially    came.       It    should  be   noted   that    for    a    maximum   value   of 
(PCF   =   1),    Figure   3   shows    the  free   surface   elevation   at    the    left 
boundary   never    falling   below   the  horizontal. 

An    interesting   relationship    developed    in    the  transport   as 
rotation   was    introduced    into   the   system,    resulting    in   a    rectifi- 
cation   of   the  along    channel    component    of    the   transport    (V) .      An 
examination    of   Figure    12    showed   that    the   cross -channel    transport    (U) 
periodically    oscillates    between   maximum  positive  and  negative   values 
while   the  along   channel    component    oscillated  between   zero  and  a 
maximum    in    one   direction.       It    was    easily    demonstrated   that    this 
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direction    is    strictly   a    function    of    initial    set-up.      As    the   flow 
was    initiated    in   the  positive  x-direction    the    initial   along   channel 
component    developed   in    the  negative   y-direction.      When    the   cross- 
channel    flow   reverses    direction,    along   channel    components    are 
generated    in    the  positive   y-direction.      These   components    are  just 
sufficient    to   cancel    the   components   previously    established    in    the 
opposite    direction,    resulting    in   a    rectification    in   the  along- 
channel    transport.      For    motions    of    large   time   scales,    this   process 
implies    that    important    changes    in    the  flow  are   caused  by   Coriolis 
accelerations    and  these   will   be    discussed   in    the   conclusions. 

Another    observed  feature    of  rotation   was    the   subsequent 
reduction    in    oscillation  period   for    increased  values    of   Coriolis. 
The   effect    of  rotation    in   reducing   the   spacial   scale    of   the  motion 
allows    the   system   to   oscillate   with  greater    frequency. 

It    is   apparent    from   the  previous    discussions    that    the 
effects    of    inertial    rotation    can   not    be   neglected  when    dealing   with 
large   time  scale   motion. 

3.      Combined   Effects    of   Friction   and  Coriolis 

The   remainder    of    the  test   runs    were  used   to    investigate 
how   friction   and   inertial    rotation    in    various    combinations    effect 
the   transient    motion.      Some    of   the    important    consequences    are    de- 
picted  in   Figures    13   through   17. 

The  phase   shift    difference   between   bottom  stress   and 
transport    existed  throughout   all    values    of    (PCF)    and    (PCD),   as 
evidenced   in   Figures    13   and   14 .      This  phase    lag    occurred    in   varying 
degrees   but   was   a    definite   characteristic   of    time   dependent   motion, 
vanishing    only   when    the   motion   became    critically    damped. 
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The   rectification    of   the  along    channel    flow  was    also  a  pre- 
dominate  characteristic    of    the  flow   where   inert ial   rotation   was 
considered.      Figures    13   and   14   again    show   this    relationship,    as 
modified  by    friction,    and    Figure    15    depicts    the   combined  effect    on 
the   free   surface.      It    can   be   seen    that    for    an   extreme   case    (PCD  =    1, 
PCF   =   1)    the   free   surface    only   very    slowly   slumps    to  a    zero  potential 
because  friction   retards    the   cross-channel    flow    long   enough   for    ro- 
tation  to  reverse    its    direction   before   the  free   surface   reaches    the 
hor  izontal . 

Figure   16    indicates    that    with  rotation    the   energy    level    is 
damped  with   relatively    fewer    oscillations.       It   must   be   realized 
however    that    as    (PCF)    increases,    the  model    is    describing   motion    on 
an    increasing    time   scale. 

Figure   17   represents    the  path   that   would  be    taken  by   fluid 
particles    (or    suspended  matter)    at    three    levels    in    the   column    over 
five  periods    of   oscillation.       In    addition    to   the  rectification  pre- 
viously   noted,    there  appears    to  be  a   net    drift    of   the   motion    in    the 
initial    cross-channel    transport    direction.      A   closer    examination    of 
Figures    13   and   14   also  show   this    drift   and    is    a    combined  result    of 
rotation   and  friction.      Since   the  motion    is    being    continuously 
damped  by    friction,    the   cross-channel    transport    can    not   return    the 
fluid  to    its    initial  position   after    each    oscillation   and  a    net 
transport    in    the    initial    direction   of   motion   results.      The   net 
drift    is    then   strictly   a    function    of   initial   set-up    conditions    of 
the   free   surface.      The  velocity   profiles    depicted    in   Figures    18 
through    20    show   some   of   the   fine   details    of   the   flow   with   the 
combined    influences    of   friction   and    inert ial   rotation. 
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Figure  18.  Effects  of  Friction  and  Inert ial 
Rotation  on  Mid-channel  Velocity 
(PCF    =    1,    PCD   =    1) . 
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Figure   19.    Effects    of   Friction   and   Inertial   Rotation 
..    on  Mid-channel   Velocity  (PCF   =    1,    PCD  =   1, 
t    =    3T/4). 
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Figure  20.    Effects    of   Friction   and   Inert ial   Rotation    on 

Mid-channel   Velocity    (PCF   =   1,    PCD   =    1,    t    =   T) 
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The  controlling  influence  on  the  stability  of  the  water 
column  continues  to  be  the  friction,  however  the  interactions  of 
the  fluid  are  now  more  complicated  due  to  the  introduction  of 
inertial  rotation.   Under  the  influence  of  increasing  frictional 
dominance  towards  the  bottom,  the  left  rotation  induced  in  the  flow 
by  Coriolis  is  increasingly  damped,  resulting  in  the  familiar  Ekman 
spiral.   During  periods  of  the  motion  previously  described  as 
smooth,  the  spiral  remains  fairly  constant  as  the  column  rotates 
to  the  left.   However,  during  the  transition  periods  (Figure  19) 
when  frictional  influence  is  reduced  due  to  decreasing  flow  velocity, 
the  spiral  unwinds  from  the  bottom  and  var ing  amounts  of  twisting  of 
the  column  results.   This  is  again  evidence  of  the  transitional 
nature  of  the  fluid  during  these  periods.   It  should  be  noted  that 
the  velocity  scale  of  Figure  19  is  an  order  of  magnitude  smaller 
than  that  of  Figure  18,  indicating  how  small  the  velocities  of  the 
column  are  during  transition.   As  the  flow  once  again  becomes  uni- 
form, the  spiral  is  quickly  re-established  under  the  stabilizing 
influence  of  bottom  friction,  as  depicted  in  Figure  20. 

It  is  interesting  to  note  from  Figure  18  that  during  the 
transition  periods  the  along  channel  component  of  velocity  does  in 
fact  change  direction  for  brief  periods.   Because  these  periods  are 
of  such  short  duration,  this  velocity  reversal  along  the  channel  does 
not  contribute  to  any  significant  transport,  in  this  direction. 

C.   DISCUSSION  OF  THE  GEOPHYSICAL  SIMULATIONS 

All  of  the  previously  discussed  processes  characteristic  of 
time  dependent  motion,  were  observed  during  runs  26  through  30. 
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However,  since  the  magnitude  of  these  processes  are  Coriolis  and 
frictional  dependent,  they  occurred  on  a  much  smaller  scale  than 
observed  for  the  test  cases. 

The  model  generated  the  time  dependent  variables  produced  by 
the  seiche  on  a  time  scale  equivalent  to  the  resonance  period  of 
the  basin.   For  these  geophysical  simulations,  the  periods  were 
all  too  small  to  allow  inertial  rotation  to  have  any  significant 
effect  on  the  motion.   However,  for  shallow  seas  of  large  hori- 
zontal extent,  values  of  (PCF)  approaching  .25  are  entirely  pos- 
sible and  test  runs  1  through  4  showed  that  inertial  effects  would 
be  significant  and  must  be  considered.   Additionally,  motion  whose 
time  scale  is  driven  by  external  forces,  independent  of  the  re- 
sonance basin,  are  likely  to  be  significantly  altered  by  rotational 
dependence.   For  example  the  long  period  waves  propagating  up  and 
down  the  continental  shelf  probably  experience  all  the  rotational 
characteristics  associated  with  the  time  dependent  motion  described 
in  the  previous  section.   The  observed  phenomenon  of  transport 
rectification  might  occur  in  these  systems  and  suggests  a  mechanism 
for  modifying  bed  load  transfer  across  the  shelf.   This  process  is 
a  function  of  initial  conditions  and  appears  to  require  a  period  of 
static  equilibrium  just  prior  to  initial  motion.   A  region  where  the 
seasonal  migration  of  atmospheric  pressure  systems  are  predominantly 
in  one  direction,  might  satisfy  these  requirements,  with  a  resulting 
net  rectification  in  one  direction.   Although  this  process  is  highly 
speculative,  it  certainly  appears  possible  and  a  modification  to 
the  model  to  facilitate  its  study  is  suggested  in  Section  V. 

The  frictional  effects  observed  during  runs  26  through  30  were 
also  minimal.   It  should  be  noted  that  Ekman  depths  of  greater  than 
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20   meters   are  possible  and   thus    these    test   runs    represent    a    minimum 
frictional    dependence. 

The   exception    to   these  minimal    effects    was    the   observed  phase 
lag  between    transport   and  bottom   stress.      This    lag   was   pronounced 
in  all    cases   and   for    the  simulation   of  South   San   Francisco   Bay, 
for    example,    averaged  40   degrees.      This    indicates    that    even   small 
amounts    of   friction   are    important   and   re-emphasizes    the   need  for    a 
mathematical    representation    of   time   dependent    transients    that    allows 
independent    calculation   of  variables    such  as   bottom   stress   and 
transport . 
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V.   CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  WORK 

The  model  demonstrated  that  eq.  (22)  was  a  suitable  mathematical 
representation  of  the  time  dependent  motion  described  by  the  governing 
equations  and  the  initial/boundary  conditions.   In  this  formulation 
the  time  dependent  variables  were  uniquely  determined  by  the  local 
time  histories  of  the  surface  at  each  point  and  this  allowed  inde- 
pendent representation  of  the  various  transients.   The  model  further 
demonstrated  the  feasibility  of  developing  a  totally  implicit  finite- 
difference  scheme  for  approximating  eqs .  (3)  and  (22)  with  accept- 
able accuracy.   With  the  insight  gained  from  the  results  it  should 
be  possible  to  improve  the  accuracy  of  the  system  by  removing  some 
of  the  important  restrictions  imposed  by  the  Ekman  dynamics  and 
the  simplified  model. 

The  physical  processes  observed  with  time  dependent  motion 
showed  the  relative  importance  of  friction  and  inertial  rotation 
in  determining  the  nature  of  the  variables.   It  was  apparent  that 
for  the  larger  time  scales  (PCF  on  the  order  of  .25  or  larger), 
Coriolis  accelerations  are  important  to  the  flow  and  must  be  con- 
sidered.  However,  shorter  time  scales  will  probably  be  adequately 
represented  neglecting  rotation  with  little  accuracy  loss.   During 
the  data  analysis  it  appeared  that  friction,  even  in  small  amounts, 
exerted  considerable  influence  on  the  characteristics  of  the  vari- 
ables and  their  relationships,  and  should  therefore  not  be  neglected 
in  representing  time  dependent  motion. 
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Since  frictional  effects  are  predominant  in  the  system,  the 
simplified  assumption  of  a  constant  coefficient  of  vertical  friction 
throughout  the  layer,  is  probably  one  of  the  most  limiting  factors 
on  the  model's  accuracy.   A  technique  for  improving  the  repre- 
sentation of  the  influence  of  friction  might  be  the  assumption  of 
constant  stress  layer  (Prandl  layer)  underlying  the  Ekman  layer. 
The  solutions  in  these  layers  would  be  matched  at  their  interface; 
that  is,  at  some  height  above  the  bottom  determined  as  a  function 
of  the  amount  of  friction  present  and  the  roughness  of  the  bottom 
boundary.   This  method  has  been  used  in  air -sea  boundary  layer 
models  with  considerable  improvement  in  accuracy,  and  merits  con- 
sideration at  the  ocean  boundary  in  this  scheme. 

The  assumption  of  constant  depth  does  not  account  for  modifi- 
cations of  the  flow  caused  by  changes  in  the  potential  vorticity 
induced  by  the  shrinking  and  stretching  of  the  water  column  over 
varying  topography.   This  restriction  is  probably  most  serious  in 
shallow  regions  where  time  dependent  flow  is  significant  and  changes 
in  the  bottom  extreme,  for  example  the  flow  in  the  vicinity  of  the 
Monterey  Submarine  Canyon. 

Certain  large  scale  motions  merit  representation  and  investi- 
gation using  this  formulation.   One  of  these  is  the  study  of 
seiches,  especially  the  trapped  modes,  on  the  continental  shelf 
since  their  buildup  can  cause  storm  surge  (Pierson  and  Neuman, 
1966).   Also,  as  previously  mentioned,  the  large  scale  oscillations 
on  the  continental  shelf  might  make  important  contributions  to 
sediment  transfer  across  the  shelf.   Therefore  it  might  prove  use- 
ful to  remove  the  horizontal  boundary  conditions  from  the  model 
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and  drive  the  flow  by  an  externally  specified  time  dependent 
forcing  term.   This  would  indicate  how  effectively  the  motion  was 
coupled  to  these  terms  and  show  whether  the  previously  noted  trans- 
port rectification  is  of  sufficient  magnitude  to  be  significant. 

An  examination  of  Tables  I  and  II  indicates  how  costly  this 
type  of  numerical  representation  is  in  terms  of  computer  time.   It 
was  not  determined  how  much  loss  in  resolution  would  occur  by  in- 
creasing the  spacial  grid  size  and  reducing  the  length  of  the  local 
time  history  considered.   It  is  clear  that  considerable  reduction 
in  computing  time  would  result  and  that  such  considerations  are 
worth  investigating  before  a  similar  model  is  used  for  extensive 
production  runs. 
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APPENDIX  A 
DETAILS  OF  THE  FINITE- DIFFERENCE  SCHEME 


Consider  first  the  calculation  of  the  components  of  the  vector 


D.,  recalling  that 
3 


D.  =  GiQ*)    ,  +  o"  -  ao"  ,  +  AtA 


(Al) 


where 


a  = 


At 
4Ax 


0    0    C2/At 
0    0    D2/At 
-10    0 


(A2) 


and 


on  = 


A  = 


D.  = 
D 


/d? 


(A3) 


(A4) 


(A5) 


Expanding  eq.  (Al)  using  (A2)  through  (A4)  yields 
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At 

3         4Ax 


0         O         C2/At 
0         0         D2/At 
-10  0 


3 

vn 

3 


At 

4ax 


0 

0 

C2/At 

0 

0 

D2/At 

-1 

0 

0 

«"-! 


n 


vj-i 


(A6) 


or 


D.    = 
3 


(C2/4Ax)(?^  +  1    -    5^)    +    U%   CI 
(D2/4Ax)(§^  +  1,-   Sj.x)    +  v"   + 

(-At/4Ax)(U^+1-    U^)    *   I1] 


(A7) 


Now   consider    the    calculation    of    the  matrix 


E.    = 
3 


eD      eD      eD 
11       12       13 


eD       e3      eD 
21       22       23 


eJ      e^      eD 
31      32      33 


,    J  ^   1 


(A8) 


Recall    eq.    (49)    and   express    it    in    inverse   form  as 


Ej  =   (I  +  aEj_1)*1  0£,   j  *  1 


(A9) 


where  I  is  the  (3X3)  identity  matrix.   Using  the  appropriate 

values  for  I,  a,  and  E.  ,,  it  is  easily  showed  that 

3-1  * 
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(I  +  aEj._1)   = 


Ate 


3-1 

11 


4ax 


Ate 


12 


4ax 


Ate 


l  - 


j-ll 

13 


4ax 


=  y 


j-i 


(A10) 


Therefore, 


(I    +  «Ej_1) 


-1 


DETCy-5-   ) 


AOJCy3"1) 


(All) 


Where    DET(y         )    is    the   determinant    of    the   matrix    (y^       )    and 
ADJ(y         )    is    its    adjoint.      The   computation   of   eq.     (All)    are   ex- 
tremely   lengthy   and  are   not    included.      The   resulting  product    of 
applying    eq.     (All)    to   eq.     (A9)    is    as    follows: 


Ej    "   DET 


j-l 


C2e 


33 


4Ax 


D2e 


33 


4ax 


A 


+C2e 


j-l' 
31 


4ax 


0 


0 


C2e 


3-1 


C2 

- 

"13 

At 

4 

Ax 

D2 

- 

D2 

•S1 

At 

4 

Ax 

cz4] 

■i  \ 

\ 

4ax 


j  *  i 


(A12) 


wher  e , 

DET.    . 
3-1 


=    1    - 


Z^Ate^-rae^V 


C2At  j-1    j-1       j-1    j-1 

(4AX)2(G33    Gl1    "G31    Gl3    ' 


(A13) 


This    completes    the   calculations    of    the   components    of   E . . 
Finally   consider    the   calculations    of   the   components    of    the   vector 
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F.  = 
D 


(A14) 


Recall  eq.  (50)  and  express  in  its  inverse  form  as 
F.  =  (I  +  «E.  ,  )_1(D.-  Q!F.  ,  )  ,  j  >    l 


(A15) 


Again  the  computations  of  eq.  (All)  are  involved  and  are  not 
included.  The  results  of  applying  eq.  (All)  to  (A15)  are  as 
follows : 


F  .  =  — — 
j    DET , 


j-l 


1  -  Ate 


u1]  Q1"  (c2dillQ3 

"5  Ax       /  V~4"ax 


lei"1    +   D2At(e^'1e^1-e^;1e?71)  ^ 


31 


4ax 


'31      13 
(4Ax) 


'e^1     ■   C2At(e^1e^1-e^1e^1) 


D2e 


"31 

ZfAx 

j-l 
33 


4ax 


Q3 


(4Ax)' 


(A16) 


,    j*l 


Ate 


j-l 
11 


k  Ax 


Ql 


1    +   C2e 


j-l 
31 


4ax 


Q3 


where, 


C2f 


j-l 


Ql    =    d^    - 


Q2    =   d^    - 


Q3  =  &> 


4ax 


d24' 

■l 

4  Ax 
Atf^" 

•  l 

4ax 

(A17) 
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Note    that    the  middle   terra   of    the   vector    F.    is    the    entire   quantity 
in   the  brackets. 

The    implicit    scheme   for    computing   U,    V,    and   §    can    now   be 
constructed.      Writing    eq.    (46)    as 

0n+!    =    E.    ,9n+1    +   P.    _  (A18) 

j-1  j-1    j  j-1 

and  substituting   the   quantities   previously    derived   for    E.    .    and 
F  .    .    yields : 

«£   -   eJ-V+1    ♦   eJ"1  5^   ♦  i{*  (A19) 

«T-l  -  4>T  +  e^  s-+1  +  *2_1  (A20) 

-n+l   eJ-l„n+l  +  J-1  f"+l  +  fJ"l  fA21x 

5j-l   e31  J     e33  &j     f3  (A  X) 

Equations  (A7) ,  (A12),  (A13),  (Al6),  (A17) ,  and  (A19)  through 
(A21)  represent  the  expressions  used  to  implicitly  calculate  values 

for  9.  at  each  grid  point  in  the  manner  described  in  the  text. 

3 
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APPENDIX   B 

■X-  ¥r 

DETAILS   OP   REPRESENTING    Cl     ,    C2     ,    VELOCITY   AND    BOTTOM    STRESS 

Consider    the  representation    of   Cl      by    first    noting    that    the 
limits    of    integration   have   the   following   meaning: 
t   =   0   represents    the    initial    condition   time. 
t    =   t   represents    the  present   time. 

t    =   t   +    At    represents    the   time  at    which   the   value   of   the 
various    variables    are  to  be    calculated. 

The    integration    of    eq.    (48)    was   performed   numerically   by 
dividing    the   total    time    into   small    intervals    (At)    and    evaluating 

p*(tM    as    a    constant   at    the  mid-point    of    each    interval.      This 

* 
allowed  Cl      to  be   numerically   represented  as    follows: 


l  h_,  _2  h_,  [  bxj 


2 
-kp    At    -if At 
e  e  -    1 


M=l      p      N=l  \        /L 

r(I-l-N)At         .^2  -ifft-tM 

X      \  -    e   *P    (t    X    }e    lt(t    X    }dt'  (Bl) 

J(I-N)At 

With   the   further    assumption   that    the    effect    of    inertial    ro- 
tation was    constant    in    each    interval,    e  was   approximated 

*    *v.         -,*         •  -if(I-^-N)At 

at    the  mid-point   as    e         x  '         and   treated   as    a    constant. 

Analytical    integration    of    the  remaining    integral    resulted    in 

eq.    (51).      To   solve   for    Cl    and  Dl    the   following   relation   was    used, 

and    the  real    and    imaginary    terms    collected. 

elft    =   cos (ft)    -    i   sin(ft)  (B2) 
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This    allowed  CI    and   Dl    to   be   represented  as 

2 


MAXM  1-1 

CI1    =   Al   S  -f  S         A3 

M=l      p      N=l 


-    cos    (f(i-^-N)At} 


MAXM  1-1 

Dl1    =   Al  S  ^§  £        A3 

M=l      p      N=l 


;"kp    Atcos{f(I+%-N)At} 


a4 


-kp    At    . 


(B3) 


sin{f(I+^-N)At} 


-   sin   {f(I-^-N)At}      A4 


(B4) 


where 


Al 


-2o. 
h 


A2    = 


kp' 


(1    " 


-kp    At) 


A3   = 


)x 


N-^ 


<£> 


a4  = 


-kp    (I-l-N)At 


Using    the   same   assumption    regarding    constant    inertial   rotation 
in    each   time   interval,    C2      readily   reduced   to   eq.    (52).      Applying 
eq.    (B2)    to   eq.     (52)    yields: 

MAXM 


C2    =   Al    cos (^f At)    E 

M=l      p 


A2 
2 


(B5) 


MAXM 


A2 


D2    =    -    Al    sin(^fAt)    E 

M=l      p 


(B6) 
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Consider  now  the  scheme  formulated  to  recover  the  velocity 
profile.   Since  velocity  was  calculated  at  time  t  +  At,  eq.  (26) 
was  represented  by: 


w 


_  2g  ^^  (-l)mcos(pz)(ct  e-kp2(t  +  At-t')e-if-(t  +  At-f)  6g 
h  M=l        P       Oo  '  bx 


dt 


X't+At  e-kp2(t  +  At-f)e-if(t  +  At-t»)  b£ 


s. 


bx 


( t ' ) dt ' 


(B7) 


Applying  the  same  assumptions  and  relations  to  eq.  (B7)  as 

were  applied  to  eq.  (Bl)  yields: 

MAXM      /    1-1     r   .2  -, 

u1  =  -  Al  S     ^   A2  S    A3   e"1^  At  cos{  f (1+%-N) At} A4 
L  M=l   P   I     N=l 


+  cos(^fAt)(A2)(A6) 


MAXM   -  /    1-1 
v.1  =  -  Al  S    ^2  ■  A2  S    A3 
L         M=l   P  C    N=l 


-e"kp  At  sin{f(I+^-N)At} 


A4 


-  sin(^fAt)(A2)(A6) 


) 


where 


M 


A5  =  (-1)   cos(pz) 


A6  =  (I1) 
v  Ox' 


i-k 


(b8) 


(B9) 


The  term  (A5)  represents  the  depth  dependence  and  (A6)  the 
slope  in  the  most  recent  time  interval.   (A6)  was  calculated  and 
stored  by  the  model  after  U,  V,  and  §  were  computed. 

The  final  quantity  needing  representation  is  bottom  stress  (T) 
Recalling  eq.  (54)  and  noting  that  the  integral  is  identical  to 
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that    of   eq.     (26),    bottom    stress    is    easily   represented   as: 


<Tx»L   = 


MAXM 
kAl  2 
M 


AXM   (        1-1  [      .2 

\k2  S        A3    U'^5        cos{f(I+3g-N)At} 
=  1     ^         N=l  L 


a4 


■} 


cos(%fAt)(A2)(A6) 


MAXM  /         1-1  w,2A+ 

(T    )       =    -    kAl  S  JA2  E         A3      -e'^5         sin{f  ( I+^-N) At} 

y  L  M=l  L 


a4 


(BIO) 


-   sin(?§fAt)(A2)(A6) 


(Bll) 


Where    (T    )    and    (T    )    are   cross-channel   and   along-channel   components 
x  y 

of   stress   respectively. 

During    the   calculation    of   CI,    Dl ,    C2 ,    D2 ,    the  model    stores 
those   common    quantities    needed    for    velocity    and   stress    calculations 
as    defined  by    eqs .    (B3)    through    ( B6 ) ,    and    (B8)    through  (Bll ) .      The 
program   recalls    these    quantities    and  applies    them   to   eqs.    (B8) 
through    (Bll)    whenever    it    is    required  to    calculate   velocity   and 
bottom  stress. 
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APPENDIX  C 
COMPUTER  PROGRAM  DESCRIPTION 
The  computer  program  was  written  in  FORTRAN  IV  and  was  used 
with  the  IBM  36O/67  computer  system  at  the  W.R.  Church  Computer 
Center,  Naval  Postgraduate  School.   The  program  has  been  generalized 
such  that  basin  size,  run  time,  initial  free  surface  elevation,  and 
number  of  modes  considered  in  the  solution  can  be  varied  in  the  data 
statements.   Additionally,  the  model  allows  variable  grid  spacing, 
requiring  only  changes  in  the  imput  parameters  and  DIMENSION  state- 
ments.  The  influence  of  inertial  rotation  and  vertical  friction  are 
introduced  into  the  model  by  the  two  non  dimensional  ratios  dis- 
cussed in  the  text. 

FORTRAN  IV  symbols  used  for  the  major  model  parameters  are 
defined  below. 
A       Model  time  (real  time  normalized  to  the  free  oscillating 

seiche  period) 
DH      Spacial  step  in  depth  (fixed  at  one  meter) 
DX      Spacial  step  in  across -channel  direction 
DT      Time  step 
ETA     Free  surface  elevation 
ETAX    Sea  surface  slope 
ETAXS    Array  for  storing  the  time  history  of  the  sea  surface 

slopes 
F       Coriolis  parameter 
FREQ    Frequency  of  the  free  oscillating  seiche 
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G  Acceleration    of   gravity 

H  Depth    of   water 

HI  Free    surface   elevation   at    left   boundary 

MAXH  Integer    representation   of  water    depth    (must    be   even) 

MAXL  Maximum  number    of    cross-channel   grid  points    (must    be    odd) 

MAXP  Interval    in    time   step    when    velocity,    energy,    and   bottom 

stress    is    calculated 

MAXT  Maximum   number    of   time   steps 

PCD  Ratio   of   Ekman    depth   to   depth    of    the  water 

PCF  Ratio   of   Coriolis   parameter    to    depth   of   the  water 

PCT  Number    of   time   steps   per    free  wave  period 

PE  Potential    energy    of   the   seiche 

PERD  Period  of  a    free   oscillating   seiche 

R  Coefficient    of   vertical    friction 

RHO  Density   of   the   water 

SE  Kinetic    energy   associated  with    the   y-component    of    velocity 

SIGMA  Angular    frequency    of   a    free   oscillating   seiche 

STRESC  y-component    of   bottom   stress 

STRESR  x-component    of   bottom  stress 

j  Real   time 

TE  Kinetic    energy   associated  with    the  x-component    of    velocity 

TOT  Total    energy    balance   of    the   seiche 

UT  x-component    of   volume   transport 

VELC  y-component    of   velocity 

VELR  x-component    of    velocity 

VT  y-component    of   volume   transport 
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The  procedure  for  initialization  of  the  program  is  as  follows: 

1)  Modify  the  initial  DATA  statements  to  correspond  to  the 
parameters  of  the  motion  being  represented. 

2)  Change  DIMENSION  statements  accordingly. 

3)  Insure  that  appropriate  storage  and  computer  time  require- 
ments are  satisfied.   The  program  required  a  maximum  of  100K  BYTES 
storage  on  the  IBM  360/67  system  and  variable  CPU  usage  time  as 
depicted  in  Tables  I  and  II. 

4)  A  minor  modification  to  the  program  is  necessary  when 
frictionless  motions  are  being  represented.   If  (PCD  =  0)  or 

(PCF  =  0),  remove  the  following  three  expressions  from  the  program 

K  =  (1-1)  -  10.0/CC 

IF(K.LE.0)K  =  1 

DO    40    J    =  K,     IM1 
Replace   these    expressions    by 

DO   40   J   =    1,    IM1 

5)  Because  of  the  formulation,  the  program  can  not  compute 
velocity  at  I  =  1;  therefore  a  value  of  MAXP  =  1  should  not  be 
used.   To  recover  the  velocity  profile  at  each  time  step  (I  j-    1), 
initialize  the  DATA  with  MAXP  =  2  and  insert  the  statement 
(MAXP  =  1)  immediately  before  statement  number  300. 

The  program,  as  appearing  in  the  next  section,  is  initialized 
in  the  MKS  system  and  changes  in  the  imput  should  be  made  ac- 
cordingly.  The  major  computational  sections  of  the  program  are 
prefaced  by  appropriate  comment  statements,  and  a  reader  with 
FORTRAN  IV  experience  should  have  little  difficulty  following  the 
program. 
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COMPUTER    PROGRAM 


DIMENSION    ETAK11)  .STA2(  11)  .ETAX( 11 )  .ETAXKlDi 
lETAXSd  1,201  )t  UT1(11),UT2(  11  )  ,VTH  11 ), VELR ( 11, 31)  t 
2VELC(11 ,31)  ,FUNR(11 ,2  5 ) ,FUNC(11  ,25 ) , P ( 25 ) , P2 ( 25 ) , 
3P3(25).P4(25)tPK(25).3(31).C(31),TIN(ll),PIN(ll), 
4P0T(11) ,E11 (10  ),i21 (10 ) ,E31(1C ),cl3(10), E23( 10), 
5E33(  10)  ,D1(  10)  ,D2( 10)  ,D3(10)  ,F1 ( 10)  ,F2 (10  )  .F3 (10  )  . 
6DET(10),CR1(  11),CC1(11),STRESR(11)  , STRE SC (  11) 

C  INPUT    INITIAL    DATA 

C 

DATA    MAXT,MAXL , MAXH , MAXP, MAXM/ 200, 11 , 30,  5 , 25/ 

DATA    WIDTH, HI , PC T, DH , RH 0/5000. .1 . #40. tl • •1027./ 

CAT A    PIE, G,PCF, PCD/ 3. 1415927, 9. 8, .5,0.0/ 
C 
C  ECHO   CHECK    DATA 

WRITE(6T310) 

WRITE (6, 320) MA  XT ,MAXL, VAXH, MAXP ,M AXM 

WRITE (6, 330 ) WI DTH, HI ,DH •RHO. PI E ,G .PCT.PCF.FCD 
C 

DX=WIDTH/FL0AT(MAXL-1 ) 

H=FLOAT(MAXH) 

PERD=(2.0*WIDTH)/SQRT(G*H) 

FREQ=1.0/P£RD 

SIGMA=2.0*PIE*FREQ 

F=PCF*SIGMA 

R=(PCD*H)**2*F/( 2.0*PIE**2) 

CT  =  PERD/PCT 

Q1=DT/(4.0*DX) 

Q2=4.0*DX 

C3=-2.C*G/H 

WRITE(6,350)H,WIDTH,DT,DX,PERD,SIGMA,F,R 
C 

C  COMPUTE    AND    STORE    INITIAL    CONDITIONS    ON    ETA ,E TAX, UT, VT 

C 

DO    10    L=1,MAXL 

X=(L-1)*DX 

ETA1(L)=HI*C0S(PIE*X/WIDTH) 

ETAX(L)  =  (-1.C*PI5*HI/WIDTH)*SIN(PIE*X/WIDTH) 

UT1(L)=0.0 

VTKL)  =0.0 
10    CONTINUE 

A  =  0.0 

ETAX(MAXL)=0.0 
C 

C  PRINT    INITIAL    CONDITIONS 

C 

WRITE(6,360)A 

WRIT £(6,370)  (ETAKL  ),L  =  1,MAXL) 

WRITE  (  6,3  80)  (ETA  X(L)  ,  L  =  l  ,  MAX  L) 

WRIT£(6,400)(UT1(L  ),L=1,MAXL) 

WRITE (6, 40 5) (VT1 (L) ,L=1,MAXL) 
C 
C  CALCULATE   AND    STORE    MODAL    CONSTANTS 


C 


DO    20    M=1,MAXM 

P(M)=(2.0*FL0AT(M)-1.0)*PIE/(2.0*H) 

P2(M)=P(M)**2 

PK(M)=P2(M)*R 

P4(M)=EXP(-PK(M)*DT) 

C5=PK(M) 

C4=C5*DT 

IF(C4.LT.1.E-4)G0    TO    15 

P3(M)=(1.0-EXP(-C4) )/C5 

GO    TO    20 
15    P3(M)=(DT*(1.0-(DT*C5/2.0)*( 1 .0-C5*DT/3.0 ) )) 
20    CONTINUE 
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c 

C  CALCULATE    CR2    ANC    CC2 

C 

FDT=p*DT 

C1=C0S(  .  5*FDT) 

C2=SIN( .5*FDT1 

CR2=0.C 

DO    30    M=1,MAXM 

CR2=CR2+P3(M)/P2(M  ) 

30  CONTINUE 
CC2=-C3*C2*CR2 
CR2=C3*C1*CR2 

C 

C  BEGIN    TIME    LOOP 

C 

DO    300    I=1,MAXT 

T=I*DT 

A=T/PERD 

IFd.NE.DGO  TO    35 

DO    31    L=1,MAXL 

CR1(L)=0.0 

CC1(L)=0.0 

31  CONTINUE 
C 

C  SET    LEFT    BCUNDRY    VALUES    FOR    MATRICES 

C 

E11(1)=0.0 

E21(l)=0.0 

E31(l)=0.0 

E13(l)=0.0 

E23( 1)=0.0 

E33(l)=1.0/.951 

Fl(l)=C.O 

F2(  1)=0.0 

F3(l)=0.0 

DET(l) =1.0 

GO    TO    6  5 
C 

C      (I.NE.l)  CALCULATE  CR1  AND  CC1 
C 

35  CO  60  L=1,MAXL 

C6=0. 0 

C7=0.0 

DO  50  M=1,MAXM 

FACR1=0.0 

FACR2=0.0 

FACC1=0.0 

FACC2=0.0 

CC=PK(M)*DT 

IM1=I-1 

K=(I-1)-10.0/CC 

IP(K.LE.0)K=1 

DO  40  J=K,IM1 

C4=ETAXS(L,J) 

C5=EXP(-CC*(I-1-J)  ) 

FACRl=FACRl+C4*C0S(FDT*U  +  .5-J)  )*C5 

FACR2=FACR2-C4*C0S( FDT*< I-. 5-J) )*C5 

FACC1=FACC1-C4*SIN(FDT*( I+.5-J ) )*C5 

FACC  2  =  FACC  2+C4*  SI N <  FDT* ( I-. 5-J ) )#C5 
40    CONTINUE 

FACR1=FACR1*P4(M) 

FACC1=FACC1*P4(M) 

C5=P3(M) 

FUNR(L,M)=FACR1*C5 

FUNC(L»M)=FACC 1*C5 

C8=1.0/P2(M) 

C6=C6+(FACR1+FACR2 )*C5*C8 

C7=C7+(FACC1+FACC2)*C5*C8 
50    CONTINUE 

CR1(L)=C3*C6 

CCKL)=C3*C7 
60    CONTINUE 
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c 

C  CALCULATE    VALUES    OF    MATRICIES 

C 

65    MM1=MAXL-1 

DO   70    L=2,MM1 

DHL)  =  (CR2/Q2)*(ETA1(L+1)-ETA1  (L-l)  )+UTl  (LI+CRKL) 

D2(L)  =  (CC2/Q2)*(ETAHL  +  1HETA1(L-1  )  )+VTl ( L ) +CC1 ( L ) 

D3(L)=-Q1*(UT1 (L+l )-UTl (L-l ) )+ETAl (L ) 

Q3=Q1/DET(L-1) 

Ell  (L)=Q3*CP 2*533  (L-D/Q2 

E13(L)=Q3*(CR2/DT-CR2*E13(L-1  )/Q2) 

E2KL  )=Q3*CC2*£33(L-1)/Q2 

E23(L)=Q3*(CC2/DT-CC2*E13(L-1)/Q2) 

E31(L)=-03M  1.0+CR2*E31(L-1)/Q2) 

E33(L  )=Q3*CR2*cll(L-l)  /Q2 

DET(L)=1.0-(1.0/Q2 )* ( DT*E1 3 ( L )-CR2*E31 (L ) )+(CR2*DT/Q2 
1**2)*(E33(L)',E11(L  )-c31(L)*E13(D) 

Q4=D1(L)-CR2*F3(L-1 )/0  2 

05=D2(L)-CC2*F3(L-1 )/Q2 

06=D3(L )+DT*FKL-l)/Q2 

Fl (L)=( (l.C-DT*E13 (L-1)/Q2)*Q4-(CR2*E33(L-1)/Q2)*Q6)/ 
lDET(L-l) 

C4=(-CC2*E31(L-1  )/02  +  CC2*DT*(E31(L-l)*E13(L-D- 
1E33(L-1)*E11(L-1 ) ) /Q2**2)*Q4 

C5=(1.0+(CR2*E31(L-1)-0T*E13(L-1) ) /02-CR2*DT*( E31 ( L-l ) 
1*E13CL-1)-E33(L-1)*E11(L-1)  )/Q2**2)*Q5 

C6=(CC2*E33(L-1)  /Q2)*Q6 

F2(L)=(C4+C5-C6)/D2T(L-1) 

F3(L)=( (DT*cll (L-l )/Q2)*Q4+(1.0+CR2*E31(L-l)/Q2)*Q6)/ 
lDET(L-l) 
70    CONTINUE 
C 

C  COMPUTE    UT2,    VT1,    ETA2 

C 

ETA2(MM1)=F3(MM1 )/(1.0-E33 (MM1 J/.951) 

ETA2(MAXL)=ETA2(MM1)/.951 

UT2(MM1  >=E13  (MM1  )*ETA2  (MAXL  ) +FKMM1 ) 

VTKMM1)  =E23(MM1  )*ETA2  (MAXD+F2  (MM1) 

UT2(MAXL )=0.0 

VT1(MAXL)=0.0 

MM3=MAXL-3 

DO    80    L=1.MM3 

K=^AXL-L 

C4=UT2(M) 

C5=ETA2(M) 

ETA2(M-1 )=E31(M-1)*C4+E33(M-1 )*C5+F3(M-1 ) 

UT2(M-1  )=£11  (M-l  )*C4+E13(M-1)*C5+F1(M-1) 

VT1(M-1)=E21(M-1 )*C4+E23(M-1)*C5+F2(M-1 ) 
80    CONTINUE 

ETA2U  )=£TA2(2)/.951 

UT2(1)=0.0 

VT1(1)=0.0 
C 
C      COMPUTE  ETAX1 


C 


DO    85    L=2,MM1 
ETAX1(L)=(ETA2(L+1)-ETA2(L-1) )/(2.0*DX) 


85    CONTINUE 

ETAX1(1)=0.0 

ETAX1(MAXL)=0.0 
C 

C  RESET    ETAXS    ARRAY   AND    RESET    ETA1  ,    ETAX,    AND   UT1 

C 

DO    90    L=1,MAXL 

ETAXS(L,I)=(ETAX1(L)+ETAX(L) )/2.0 

UT1(L)=UT2(L ) 

ETA1(L)=£TA2(L  ) 

ETAX(L)=ETAX1(L) 
90    CONTINUE 
C 

C  PRINT    OUT    UT1,     VT1,    ETA1,    AND    ETAX 

C 
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WRITE(6,410 )A 

WRITE  (6,420)  (UTKL)  tL=l,MAXL) 

WRITE(6,421XVT1(L  )tL-l,MAXL) 

WRITE (6, 42 5)  ( ET Al ( L  )  , L  =  l , MAXL ) 

WRITE( 6,429) ( E TA X( L ) . L=l . MA XL ) 
C 

IF(MOD(I ,MAXP) .EQ.O)GO    TO    100 

GO    TO    300 
C 
C  I    IS    A    MULTIPLE    OF    MAXP,    CALCULATE    VELOCITY 


C 


100    CO    115    L=1»MAXL 
C6=ETAXS(L,I  ) 
DO    110    M=1,MAXM 

FUNR(L,M)=FUNR(L»M)+C1*P3(M)*C6 
FUNCC  L,M)=FUNC(L,M)-C2*P3(M)*C6 


110 

CONTINUE 

115 

CONTINUE 

DO    140    L=1,MAXL 

MH1=MAXH+1 

DO    130     J=1,MH1 

Z=(J-1)*DH 

SI=-1.0 

C4=0.0 

C5=0.0 

DO    120    M=1TMAXM 

C6=P(M) 

C7=C0S(C6*Z)  /C6 

C4=C4+SI*FUNR(L,M)*C7 

C5=C5+SI*FUNC<  L,M)*C7 

SI=-1.0*SI 

120 

CONTINUE 

V£LR(L,J)=-C3*C4 
VELC(L, J)=-C3*C5 

130 

CONTINUE 

c 
c 
c 

140 

CONTINUE 

CALCULATE    BCTTOM    STRESS 

C6=C3*R 

DC    142    L=l ,MAXL 

C4=0.0 

C5=0.0 

DO    141     M=1,MAXM 

C4=C4+FUNR(L,M) 
C5=C5+FUNC(L,M ) 

141  CONTINUE 
STRESR(L)=C4*C6 
STRESC(L)  =  C5*C6 

142  CONTINUE 
C 

C  CALCULATE    ENERGY    BALANCE    USING    SIMPSON1 S   RULE 


C 


145    DO    165    L=ltMAXL 
DO    150    J=1,MH1 
B( J)=VELR(L,J)**2 
C( J)=VELC(L,J)**2 

150    CONTINUE 

EVSUM1=0.0 

EVSUM 2=0.0 

0DSUM1  =  0.0 

0DSUM2=0.0 

ENSUM1=B(  1)+B(  MH1) 

ENSUM2=C(1 )+C(MHl ) 

MH2=MH1-1 

DO    155    J=2,MH2,2 

EVSUM1=EVSUM1+B( J) 

EVSUM2=E VSUM2+C( J) 

155    CONTINUE 
^H2=MHl-2 
DO    160    J=3,MH2,2 
0DSUM1=0DSUM1+B(J ) 
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c 
c 
c 


c 
c 
c 


0DSUM2=0DSUM2+C( J) 
160    CONTINUE 

TIN(L)=0H/3.0*(ENSUM1+2.0*0DSUM1+4.0*EVSUM1) 

PIN(L)=DH/3.0*  (ENSUM2+2.0*0DSUM2+4.0*£VSUM2) 

P0T(L)=ETA1(L)**2 
165    CONTINUE 

EVSUM1=0.0 

EVSUM 2=0.0 

EVSUM3=0.0 

0DSUM1=0.0 

0DSUM2=0.0 

0CSUM3=0.0 

ENSUM1=P0T(1)+P0T( MAXL) 

ENSUM2=TIN( 1  )+TIN(MAXL ) 

ENSUM3  =  PIN(1 )+PIN(MAXL) 

DO    170    L=2,MM1,2 

EVSUM1=EVSUM1+P0T(L) 

EVSUM2=EVSUM2  +  TIN(  L) 

EVSUM3=EVSUM3+PIN(  L) 
170    CONTINUE 

MM2=MAXL-2 

DO    175    L=3,MM2,2 

0DSUM1=0DSUM1+P0T( L) 

0DSUM2=0DSUM2+TI N( L) 

0DSUM3=0DSUM3+PIN(L) 
175    CONTINUE 

PE=.5*RH0*G*(DX/3.  0)*(ENSUM1+    2  .  0*0DSUM1  +4  .0*EVSUM1  ) 

TE=.5*RH0*(DX/3.0)*(£NSUM2+2.0*ODSUM2+4.O*EVSUM21 

SE=.  5*RH0*  (DX/3.0)*(ENSUM3+2.0#0DSUM3+4.0*£VSUM3) 

TOT=P£+TE+SE 

PRINT    VELOCITY 

WRITE(6,430  ) 

DO    180    J=1,MH1 

JM1=J-1 

WRITE (6, 440) JM1, ( V5LR ( L, J ) , L=l , MAXL ) 
180    CONTINUE 

WRITE(6,435) 

DO    190    J=1,MH1 

JM1=J-1 

WRITE(6,440)JM1,  (VELC(L,J)  ,L  =  1,MAXL) 
190    CONTINUE 

PRINT    PE,TE,SE,TOT,STRESR,    AND    STRESC 

WRI TE< 6,45  0) PE ,TE,  SE,TOT 
WRITSC6. 45  5) ( STRESR(L )  ,L=1, MAXL) 
WRITE (6 ,460) (STRESC (L) ,L=1,MAXL  ) 


300    CONTINUE 


STOP 

310    FORM  AT ('l1 

320    FORMAT(»0» 

330    FORMAT( '0' 

350    FORMAT  CI1 

l'FGLLOWSS/aX,'  D 
2'WIDTH=',F7. 1, IX, 


•ECHO    CHECK    OF    INPUT    DATA',/) 


»  s 
.9 


360 
] 
370 

3  80 
400 
405 
410 

4  20 
4  21 
425 
429 


31X,,DX=« ,F6.1, IX 
4/,lX,«  SIGMA=«  ,E11 


11.  3) 

•SOLUTION  TO 

• nsPTH=i tF5 

•METERS'  , 


SEICHE 


FORMAT ( 
F4.1 ,1X 
FORMAT( 
FORMAT ( 
FORMAT ( 
FORM£T( 
FORMAT ( 
FORMAT( 
FORMAT ( 
FORMAT( 
FORMAT( 


PROBLEM 
,., "METERS 
1X,'DT  =  '  »F5 
•METERS  •,/,  IX,  «PERD=», 
.3,/,lX,lF=«  , Ell. 3,/, IX, • 
CONDITIONS', /,56X, 


WITH    CORIOLI  S* 


,/,  IX, 
1.1X.1 


SEC  ,/ 


55X, •  INITIAL 
•SECONDS' ) 

•ETA'/( 11E11.3) ) 

•ETAX'/(  11E11.3)  ) 
55X,'TRANSPCRT(UT)'/ (11  El 1.3) 


,55X 
t55X, 


55X 
55X 


•  TRANSPORT*  VT)  • 


NORMALIZED    T  IME= ' , F8 .2, / ) 


/(11E11 
i 


3)  ) 


R=',E11.3,// 
•TIME=I , 


,55X,' TRANSPORT (UT) »/ (11511.3) ) 
,55X, • TRANSPORT! VT)  •/( 11E11.3) ) 
,55X  v'ETAa/(llE11.3) ) 
,55X, •ETAX,/(llcll.3) ) 


88 


430  FORMAT  CO 
435  FORMATCO 
440  FORMAT* * 
450  FORMAT ( '0 
1E15.7,8X, 
455  FORMAT( «0 
460  FORMATCO 
END 


•VELOCITY! VELR) »,/, 3X, 'Z' 

t  «Z« 
2.(  11=11.3)  ) 
■t 2Xt E15«7i  8X, 
,E15.7) 
•STRSSS(R)  ■ 


5  5X 
f55X,« VELO 

»  2X.I 
,  »PE= 

TOT=' 
55X 


»/) 
,/) 


•TE=,,E15.7,8X, «SE= 


55X, 'STRESS(C) • 


11F11.3) ) 
11E11.3)) 
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